The Role of Self-Regulation on Resiliency after TBI

Set up

Import Data

Tables by Severity

Demographic

The table below provides demographic information based on severity of injury

# level naming for categorical variables
df_demo$gender <- factor(df_demo$gender,
                   levels = c(1,2,3),
                   labels = c("Male", "Female", "Nonbinary"))

df_demo$work_current <- factor(df_demo$work_current,
                   levels = c(1,0),
                   labels = c("Yes", "No"))

df_demo$severity <- factor(df_demo$severity,
                   levels = c(2,3),
                   labels = c("Moderate", "Severe"))

df_demo$mech_injury <- factor(df_demo$mech_injury,
                   levels = c(1,2,3,4,5),
                   labels = c("Fall", "MVC", "Sports", "Violence", "Pedestrian struck"))

df_demo$income <- factor(df_demo$income,
                   levels = c(1,2,3),
                   labels = c("<52K", "52K-156K", ">156K"))

df_demo$marital_status <- factor(df_demo$marital_status,
                                 levels = c(1, 2, 3, 4),
                                 labels = c("Single", "Married", "Divorced", "Widowed"))

demo_table <- df_demo %>%
  subset(., select = c(age_current, time_injury, gender, edu, race, work_current, income, house_size, marital_status, substance, severity, mech_injury)) %>%
  tbl_summary(
    missing = "no",
    by = severity,
    type = list(
      c(age_current, edu, house_size, substance, time_injury) ~ "continuous",
      c(gender, income, work_current, severity, mech_injury) ~ "categorical"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} ({p}%)"),
    label = list(
      age_current ~ "Age (years)",
      time_injury ~ "Time since TBI (years)",
      gender ~ "Gender",
      race ~ "Race/Ethnicity",
      edu ~ "Education (years)",
      work_current ~ "Employment status",
      income ~ "Annual household income",
      house_size ~ "Size household",
      marital_status ~ "Marital status",
      substance ~ "Substance use score",
      mech_injury ~ "Cause of injury",
      severity ~ "Severity of Injury"
    )
  ) %>%
  add_p(
    test = list(all_continuous() ~ "t.test", all_categorical() ~ "chisq.test"),
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  add_n()

print(demo_table)   

Characteristic	N	Moderate, N = 191	Severe, N = 261	p-value2
Age (years)	45	51 (14)	44 (14)	0.10
Time since TBI (years)	45	7 (5)	10 (8)	0.19
Gender	45			0.24
Male		7 (37%)	15 (58%)
Female		11 (58%)	11 (42%)
Nonbinary		1 (5.3%)	0 (0%)
Education (years)	45	15.47 (2.04)	15.15 (2.65)	0.65
Race/Ethnicity	45			0.20
Asian		0 (0%)	2 (7.7%)
Biracial		2 (11%)	0 (0%)
Hispanic		1 (5.3%)	3 (12%)
White		16 (84%)	21 (81%)
Employment status	45			0.77
Yes		9 (47%)	10 (38%)
No		10 (53%)	16 (62%)
Annual household income	45			0.84
<52K		6 (32%)	10 (38%)
52K-156K		9 (47%)	12 (46%)
>156K		4 (21%)	4 (15%)
Size household	45	2.00 (1.05)	2.19 (1.39)	0.60
Marital status	45			0.14
Single		5 (26%)	14 (54%)
Married		11 (58%)	8 (31%)
Divorced		3 (16%)	4 (15%)
Widowed		0 (0%)	0 (0%)
Substance use score	45	4.16 (3.62)	1.73 (1.93)	0.013
Cause of injury	45			0.14
Fall		10 (53%)	5 (19%)
MVC		4 (21%)	10 (38%)
Sports		1 (5.3%)	4 (15%)
Violence		1 (5.3%)	4 (15%)
Pedestrian struck		3 (16%)	3 (12%)
1 Mean (SD); n (%)
2 Welch Two Sample t-test; Pearson’s Chi-squared test

ACS

ACS3 (activity re-engagement scores - outcome measure) by severity of injury

acs_table <- df_demo %>%
  subset(., select = c(acsg_prev, acsg_curr, acsg_retain, acsi_prev, acsi_curr, acsi_retain, acsl_prev, acsl_curr, acsl_retain, acsf_prev, acsf_curr, acsf_retain, acss_prev, acss_curr, acss_retain, severity)) %>%
  tbl_summary(
    missing = "no",
    by = severity,
    type = list(
      c(acsg_prev, acsg_curr, acsg_retain, acsi_prev, acsi_curr, acsi_retain, acsl_prev, acsl_curr, acsl_retain, acsf_prev, acsf_curr, acsf_retain, acss_prev, acss_curr, acss_retain) ~ "continuous",
      c(severity) ~ "categorical"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} ({p}%)"),
    label = list(
      acsg_prev ~ "ACS Global Before",
      acsg_curr ~ "ACS Global Current",
      acsg_retain ~ "Global Retained (%)",
      acsi_prev ~ "ACS IADL Before",
      acsi_curr ~ "ACS IADL Current",
      acsi_retain ~ "IADL Retained (%)",
      acsl_prev ~ "ACS Leisure Before",
      acsl_curr ~ "ACS Leisure Current",
      acsl_retain ~ "Leisure Retained (%)",
      acsf_prev ~ "ACS Fitness Before",
      acsf_curr ~ "ACS Fitness Current",
      acsf_retain ~ "Fitness Retained (%)",
      acss_prev ~ "ACS Social Before",
      acss_curr ~ "ACS Social Current",
      acss_retain ~ "Social Retained (%)",
      severity ~ "Severity of Injury"
    )
  ) %>%
  add_p(
    test = list(all_continuous() ~ "t.test", all_categorical() ~ "chisq.test"),
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  add_n()

print(acs_table)   

Characteristic	N	Moderate, N = 191	Severe, N = 261	p-value2
ACS Global Before	45	72 (11)	68 (9)	0.23
ACS Global Current	45	54 (17)	52 (13)	0.66
Global Retained (%)	45	75 (19)	76 (15)	0.77
ACS IADL Before	45	22.16 (2.27)	21.04 (3.33)	0.19
ACS IADL Current	45	17.9 (4.9)	16.8 (4.5)	0.42
IADL Retained (%)	45	81 (19)	80 (18)	0.91
ACS Leisure Before	45	22.5 (5.8)	21.1 (4.5)	0.39
ACS Leisure Current	45	18.0 (6.0)	16.8 (5.0)	0.46
Leisure Retained (%)	45	82 (22)	80 (16)	0.75
ACS Fitness Before	45	13.2 (4.4)	13.0 (4.4)	0.85
ACS Fitness Current	45	8.3 (4.7)	8.4 (3.4)	0.92
Fitness Retained (%)	45	64 (32)	69 (35)	0.61
ACS Social Before	45	13.95 (1.22)	12.88 (1.58)	0.015
ACS Social Current	45	9.68 (3.08)	9.94 (2.52)	0.77
Social Retained (%)	45	69 (19)	77 (18)	0.15
1 Mean (SD)
2 Welch Two Sample t-test

Below is the ttest for the specific t and p value for the difference between ACS3 previous social score, which was significantly different.

#significant difference in ACSS between groups. t-test for t value
df_mod <- subset(df, severity == 2)
df_severe <- subset(df, severity == 3)
ttest_acsstotal <- t.test(df_mod$acss_prev, df_severe$acss_prev)
print(ttest_acsstotal)

## 
##  Welch Two Sample t-test
## 
## data:  df_mod$acss_prev and df_severe$acss_prev
## t = 2.5391, df = 42.829, p-value = 0.01483
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2185693 1.9069368
## sample estimates:
## mean of x mean of y 
##  13.94737  12.88462

FrSBe

Comparison of self-regulation scores by severity

frsbe_table <- df_demo %>%
  subset(., select = c(frsbe_exec, frsbe_disinhib, frsbe_apathy, frsbe_total, severity)) %>%
  tbl_summary(
    missing = "no",
    by = severity,
    type = list(
      c(frsbe_exec, frsbe_disinhib, frsbe_apathy, frsbe_total) ~ "continuous",
      c(severity) ~ "categorical"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} ({p}%)"),
    label = list(
      frsbe_exec ~ "Executive function",
      frsbe_disinhib ~ "Disinhibition",
      frsbe_apathy ~ "Apathy",
      frsbe_total ~ "Total FrSBe Score",
      severity ~ "Severity of Injury"
    )
  ) %>%
  add_p(
    test = list(all_continuous() ~ "t.test", all_categorical() ~ "chisq.test"),
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  add_n()

print(frsbe_table)   

Characteristic	N	Moderate, N = 191	Severe, N = 261	p-value2
Executive function	45	41 (10)	42 (11)	0.93
Disinhibition	45	32.3 (6.1)	33.0 (6.0)	0.71
Apathy	45	33 (8)	33 (9)	0.82
Total FrSBe Score	45	106 (19)	108 (21)	0.80
1 Mean (SD)
2 Welch Two Sample t-test

TBI Composite

Comparison of composite scores for TBI QOL by severity of injury. Composite scores were calculated using:

Tyner, C. E., Boulton, A. J., Sherer, M., Kisala, P. A., Glutting, J. J., & Tulsky, D. S. (2020). Development of Composite Scores for the TBI-QOL. Arch Phys Med Rehabil, 101(1), 43-53. https://doi.org/10.1016/j.apmr.2018.05.036

tbiqol_composite_table <- df_demo %>%
  subset(., select = c(phys_health_index, emo_health_index, cog_health_index, soc_health_index, glob_health_index, severity)) %>%
  tbl_summary(
    missing = "no",
    by = severity,
    type = list(
      c(phys_health_index, emo_health_index, cog_health_index, soc_health_index, glob_health_index) ~ "continuous",
      c(severity) ~ "categorical"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} ({p}%)"),
    label = list(
      phys_health_index ~ "Physical Health Index",
      emo_health_index ~ "Emotional Health Index",
      cog_health_index ~ "Cognitive Health Index",
      soc_health_index ~ "Social Health Index",
      glob_health_index ~ "Global Health Index",
      severity ~ "Severity of Injury"
    )
  ) %>%
  add_p(
    test = list(all_continuous() ~ "t.test", all_categorical() ~ "chisq.test"),
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  add_n()

print(tbiqol_composite_table)   

Characteristic	N	Moderate, N = 191	Severe, N = 261	p-value2
Physical Health Index	45	91 (14)	97 (14)	0.15
Emotional Health Index	45	97 (12)	102 (15)	0.20
Cognitive Health Index	45	93 (13)	96 (16)	0.59
Social Health Index	45	94 (12)	91 (13)	0.43
Global Health Index	45	93 (13)	96 (14)	0.42
1 Mean (SD)
2 Welch Two Sample t-test

Predictive variables

Table 2 in dissertation

This table compares only the Personal and Environmental Protective factors and self-regulation outlined in the dissertation. Note that the Cognitive Health Composite score was not used as it includes executive functioning, which in this paper is considered a self-regulatory process. Therefore, general cognitive functioning was used which assesses memory and concentration.

PPF_table <- df_demo %>%
  subset(., select = c(phys_health_index, emo_health_index, tbiqol_genconcern_tscore, bfi_extraversion, bfi_agreeable, bfi_consciousness, bfi_neuroticism, bfi_openness, income, marital_status, spstotal, frsbe_exec, frsbe_disinhib, frsbe_apathy, frsbe_total, severity)) %>%
  tbl_summary(
    missing = "no",
    by = severity,
    type = list(
      c(phys_health_index, emo_health_index, tbiqol_genconcern_tscore, bfi_extraversion, bfi_agreeable, bfi_consciousness, bfi_neuroticism, bfi_openness) ~ "continuous",
      c(severity, income) ~ "categorical"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} ({p}%)"),
    label = list(
      phys_health_index ~ "Physical Health Index",
      emo_health_index ~ "Emotional Health Index",
      tbiqol_genconcern_tscore ~ "General Cognition",
      bfi_extraversion ~ "Extraversion", 
      bfi_agreeable ~ "Agreeable", 
      bfi_consciousness ~ "Consciousness", 
      bfi_neuroticism ~ "Neuroticism", 
      bfi_openness ~ "Openness",
      income ~ "Annual household income",
      marital_status ~ "Marital status",
      spstotal ~ "Social Support",
      frsbe_exec ~ "Executive function",
      frsbe_disinhib ~ "Disinhibition",
      frsbe_apathy ~ "Apathy",
      frsbe_total ~ "Total score",
      severity ~ "Severity of Injury"
    )
  ) %>%
  add_p(
    test = list(all_continuous() ~ "t.test", all_categorical() ~ "chisq.test"),
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  add_n()

print(PPF_table)   

Characteristic	N	Moderate, N = 191	Severe, N = 261	p-value2
Physical Health Index	45	91 (14)	97 (14)	0.15
Emotional Health Index	45	97 (12)	102 (15)	0.20
General Cognition	45	36 (8)	37 (9)	0.57
Extraversion	45	7.16 (2.50)	6.96 (2.13)	0.78
Agreeable	45	7.11 (1.94)	7.27 (2.05)	0.79
Consciousness	45	8.16 (1.54)	7.73 (1.93)	0.41
Neuroticism	45	6.47 (2.20)	6.00 (2.59)	0.51
Openness	45	8.53 (2.09)	7.42 (1.88)	0.076
Annual household income	45			0.84
<52K		6 (32%)	10 (38%)
52K-156K		9 (47%)	12 (46%)
>156K		4 (21%)	4 (15%)
Marital status	45			0.14
Single		5 (26%)	14 (54%)
Married		11 (58%)	8 (31%)
Divorced		3 (16%)	4 (15%)
Widowed		0 (0%)	0 (0%)
Social Support	45	84 (10)	76 (11)	0.022
Executive function	45	41 (10)	42 (11)	0.93
Disinhibition	45	32.3 (6.1)	33.0 (6.0)	0.71
Apathy	45	33 (8)	33 (9)	0.82
Total score	45	106 (19)	108 (21)	0.80
1 Mean (SD); n (%)
2 Welch Two Sample t-test; Pearson’s Chi-squared test

#significant difference in ACSS between groups. t-test for t value
df_mod <- subset(df, severity == 2)
df_severe <- subset(df, severity == 3)
ttest_openness <- t.test(df_mod$bfi_openness, df_severe$bfi_openness)
print(ttest_openness)

## 
##  Welch Two Sample t-test
## 
## data:  df_mod$bfi_openness and df_severe$bfi_openness
## t = 1.8232, df = 36.394, p-value = 0.07649
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1234927  2.3299705
## sample estimates:
## mean of x mean of y 
##  8.526316  7.423077

Below is the t-test for SPS total, which was significantly different between severity of injury

#significant difference in SPS between groups. t-test for t value

ttest_spstotal <- t.test(df_mod$spstotal, df_severe$spstotal)
print(ttest_spstotal)

## 
##  Welch Two Sample t-test
## 
## data:  df_mod$spstotal and df_severe$spstotal
## t = 2.3726, df = 41.268, p-value = 0.0224
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   1.13607 14.11494
## sample estimates:
## mean of x mean of y 
##  83.89474  76.26923

Descriptive Statistics

Descriptive statistics for each variable of interest for the data set including mean, median, SD, and IQR, kurtosis and se

summary_stats <- describeBy(df) %>% round(2)

# Display the summary statistics using kable
knitr::kable(summary_stats)

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
record_id*	1	45	23.00	13.13	23.00	23.00	16.31	1.0	45.0	44.0	0.00	-1.28	1.96
age_current	2	45	47.24	14.55	48.00	47.19	20.76	21.0	72.0	51.0	0.14	-1.21	2.17
age_injury	3	45	38.40	14.92	35.00	37.73	17.79	18.0	66.0	48.0	0.33	-1.28	2.22
time_injury	4	45	8.91	7.08	7.00	7.96	5.93	1.0	30.0	29.0	1.21	0.94	1.06
gender	5	45	1.53	0.55	2.00	1.51	1.48	1.0	3.0	2.0	0.28	-1.16	0.08
race*	6	45	3.69	0.76	4.00	3.89	0.00	1.0	4.0	3.0	-2.43	4.99	0.11
edu	7	45	15.29	2.39	16.00	15.32	2.97	10.0	20.0	10.0	-0.22	-0.70	0.36
work_current	8	45	0.42	0.50	0.00	0.41	0.00	0.0	1.0	1.0	0.30	-1.95	0.07
hours_work	9	19	30.50	13.04	40.00	31.50	0.00	4.0	40.0	36.0	-0.74	-1.19	2.99
occ_years_pos	10	18	7.43	9.21	3.50	6.47	5.04	0.1	30.0	29.9	1.47	1.14	2.17
diff_occ	11	19	0.53	0.51	1.00	0.53	0.00	0.0	1.0	1.0	-0.10	-2.09	0.12
no_occ_stat	12	26	3.19	0.75	3.00	3.14	0.00	2.0	5.0	3.0	1.35	1.46	0.15
income	13	45	1.82	0.72	2.00	1.78	1.48	1.0	3.0	2.0	0.26	-1.08	0.11
house_size	14	45	2.11	1.25	2.00	1.92	1.48	1.0	7.0	6.0	1.65	3.50	0.19
children	15	45	0.44	0.50	0.00	0.43	0.00	0.0	1.0	1.0	0.22	-2.00	0.07
num_child	16	20	2.00	1.26	2.00	1.75	1.48	1.0	5.0	4.0	1.21	0.49	0.28
severity	17	45	2.58	0.50	3.00	2.59	0.00	2.0	3.0	1.0	-0.30	-1.95	0.07
mech_injury	18	45	2.40	1.40	2.00	2.27	1.48	1.0	5.0	4.0	0.68	-0.91	0.21
mech_injury_other*	19	6	1.17	0.41	1.00	1.17	0.00	1.0	2.0	1.0	1.36	-0.08	0.17
substance	20	45	2.76	2.99	2.00	2.30	2.97	0.0	13.0	13.0	1.31	1.70	0.45
acsg_prev	21	45	69.60	10.14	68.00	69.43	10.38	48.0	93.0	45.0	0.25	-0.58	1.51
acsg_curr	22	45	52.70	14.33	49.50	51.85	14.08	29.5	85.5	56.0	0.51	-0.48	2.14
acsg_retain	23	45	75.80	16.61	74.00	75.97	14.83	43.0	107.0	64.0	-0.06	-0.49	2.48
acsi_prev	24	45	21.51	2.95	21.00	21.62	2.97	12.0	26.0	14.0	-0.55	0.61	0.44
acsi_curr	25	45	17.24	4.66	17.00	16.95	5.93	10.5	26.0	15.5	0.42	-1.10	0.69
acsi_retain	26	45	80.36	18.31	82.00	80.89	22.24	42.0	110.0	68.0	-0.17	-1.09	2.73
acsl_prev	27	45	21.67	5.09	21.00	21.73	5.93	11.0	32.0	21.0	-0.06	-0.81	0.76
acsl_curr	28	45	17.28	5.41	16.00	17.03	5.93	8.5	28.0	19.5	0.39	-1.02	0.81
acsl_retain	29	45	80.47	18.66	78.00	80.68	17.79	44.0	122.0	78.0	-0.02	-0.57	2.78
acsf_prev	30	45	13.07	4.35	13.00	13.19	4.45	4.0	20.0	16.0	-0.27	-0.91	0.65
acsf_curr	31	45	8.33	3.97	8.50	8.15	2.97	1.0	18.0	17.0	0.31	-0.18	0.59
acsf_retain	32	45	66.51	33.40	63.00	63.00	22.24	9.0	200.0	191.0	1.58	4.05	4.98
acss_prev	33	45	13.33	1.52	14.00	13.43	1.48	9.0	17.0	8.0	-0.44	0.38	0.23
acss_curr	34	45	9.83	2.74	10.00	9.91	2.97	4.5	15.0	10.5	-0.27	-0.66	0.41
acss_retain	35	45	73.78	19.01	73.00	74.38	19.27	35.0	109.0	74.0	-0.21	-0.89	2.83
activity_card_sort_complete	36	45	2.00	0.00	2.00	2.00	0.00	2.0	2.0	0.0	NaN	NaN	0.00
spstotal	37	45	79.49	11.38	82.00	80.05	14.83	55.0	96.0	41.0	-0.36	-1.19	1.70
bfi_extraversion	38	45	7.04	2.27	7.00	7.14	2.97	2.0	10.0	8.0	-0.18	-1.09	0.34
bfi_agreeable	39	45	7.20	1.98	7.00	7.30	2.97	3.0	10.0	7.0	-0.44	-0.81	0.30
bfi_consciousness	40	45	7.91	1.77	8.00	8.08	1.48	3.0	10.0	7.0	-0.76	-0.21	0.26
bfi_neuroticism	41	45	6.20	2.42	6.00	6.22	2.97	2.0	10.0	8.0	0.00	-1.12	0.36
bfi_openness	42	45	7.89	2.03	8.00	8.08	2.97	2.0	10.0	8.0	-0.59	-0.34	0.30
frsbe_exec	43	45	41.58	10.06	43.00	41.22	10.38	24.0	63.0	39.0	0.16	-0.77	1.50
frsbe_apathy	44	45	32.84	8.17	32.00	32.49	10.38	18.0	52.0	34.0	0.36	-0.81	1.22
frsbe_disinhib	45	45	32.67	6.01	32.00	32.51	7.41	21.0	46.0	25.0	0.20	-0.80	0.90
frsbe_total	46	45	107.09	19.77	109.00	106.46	22.24	72.0	150.0	78.0	0.23	-0.81	2.95
frsbe_complete	47	45	2.00	0.00	2.00	2.00	0.00	2.0	2.0	0.0	NaN	NaN	0.00
tbiqol_part_sra_tscore	48	45	46.40	6.76	46.10	45.88	6.08	32.1	64.1	32.0	0.80	1.01	1.01
tbiqol_anger_tscore	49	45	50.95	9.98	51.50	50.81	11.42	33.1	69.9	36.8	0.11	-1.03	1.49
tbiqol_anxiety_tscore	50	45	55.48	9.27	55.30	55.71	9.49	36.1	73.0	36.9	-0.20	-0.75	1.38
tbiqol_comm_tscore	51	45	46.58	9.39	46.20	46.51	7.86	29.2	65.5	36.3	0.15	-0.84	1.40
tbiqol_depression_tscore	52	45	53.66	9.69	53.80	53.75	10.23	33.6	74.0	40.4	-0.05	-0.67	1.44
tbiqol_dyscontrol_tscore	53	45	50.68	8.25	52.20	51.07	7.71	33.2	66.8	33.6	-0.39	-0.48	1.23
tbiqol_execfunc_tscore	54	45	35.63	6.04	34.30	35.39	5.19	24.3	50.8	26.5	0.35	-0.50	0.90
tbiqol_fatigue_tscore	55	45	54.73	8.50	54.60	54.75	8.30	37.9	72.5	34.6	0.04	-0.63	1.27
tbiqol_genconcern_tscore	56	45	36.49	8.61	36.00	36.28	8.75	19.7	53.8	34.1	0.20	-0.68	1.28
tbiqol_grief_tscore	57	45	52.42	9.48	53.40	52.95	7.26	30.7	70.3	39.6	-0.57	-0.18	1.41
tbiqol_headache_tscore	58	45	49.24	9.23	49.50	48.74	13.64	38.5	67.1	28.6	0.19	-1.31	1.38
tbiqol_mobility_tscore	59	45	45.77	8.83	44.70	45.35	8.75	31.5	63.6	32.1	0.43	-0.77	1.32
tbiqol_pain_tscore	60	45	55.23	10.85	57.40	55.20	11.27	38.4	74.8	36.4	-0.23	-1.13	1.62
tbiqol_posaffect_tscore	61	45	50.30	7.38	50.00	50.19	8.15	35.4	68.9	33.5	0.18	-0.47	1.10
tbiqol_resilience_tscore	62	45	48.70	7.99	49.10	48.39	7.41	33.4	73.6	40.2	0.44	0.53	1.19
tbiqol_selfesteem_tscore	63	45	48.38	10.63	48.90	48.48	10.53	28.4	66.0	37.6	-0.05	-0.92	1.59
tbiqol_satissra_tscore	64	45	45.28	6.31	45.10	44.85	4.60	34.7	63.2	28.5	0.85	1.09	0.94
tbiqol_stigma_tscore	65	44	50.81	7.38	51.95	51.41	5.63	33.5	62.3	28.8	-0.71	-0.14	1.11
tbiqol_ue_tscore	66	45	44.34	8.56	42.50	44.27	8.30	27.9	58.1	30.2	0.36	-0.92	1.28
marital_status	67	45	1.73	0.72	2.00	1.68	1.48	1.0	3.0	2.0	0.42	-1.05	0.11
phys_health	68	45	109.96	17.26	110.50	109.89	21.35	82.6	143.3	60.7	-0.03	-1.12	2.57
phys_health_index	69	45	94.80	13.91	95.00	95.03	16.31	64.0	117.0	53.0	-0.16	-0.97	2.07
emo_health	70	45	160.10	25.61	158.00	160.56	30.84	115.1	201.3	86.2	-0.06	-1.26	3.82
emo_health_index	71	45	99.69	14.04	101.00	99.54	17.79	77.0	123.0	46.0	0.01	-1.33	2.09
cog_health	72	45	72.12	14.14	70.70	71.81	14.97	44.0	104.6	60.6	0.28	-0.57	2.11
cog_health_index	73	45	94.60	14.51	94.00	94.68	17.79	62.0	123.0	61.0	0.00	-0.68	2.16
soc_health	74	45	91.69	12.16	90.30	91.03	10.82	69.6	127.3	57.7	0.64	0.51	1.81
soc_health_index	75	45	92.71	12.36	93.00	92.97	10.38	64.0	122.0	58.0	-0.18	0.21	1.84
glob_health	76	45	381.80	45.41	380.00	381.22	54.86	303.0	468.0	165.0	0.04	-1.11	6.77
glob_health_index	77	45	94.73	13.36	95.00	94.54	14.83	71.0	120.0	49.0	0.00	-1.03	1.99

ACS3

Outcome variable: ACS3 for all scores mean(sd)

mean_sd_acs <- df_demo %>%
  subset(., select = c(acsg_prev, acsg_curr, acsg_retain, acsi_prev, acsi_curr, acsi_retain, acsl_prev, acsl_curr, acsl_retain, acsf_prev, acsf_curr, acsf_retain, acss_prev, acss_curr, acss_retain)) %>%
  tbl_summary(
    missing = "no",
    type = list(
      c(acsg_prev, acsg_curr, acsg_retain, acsi_prev, acsi_curr, acsi_retain, acsl_prev, acsl_curr, acsl_retain, acsf_prev, acsf_curr, acsf_retain, acss_prev, acss_curr, acss_retain) ~ "continuous"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})"),
    label = list(
      acsg_prev ~ "ACS Global Before",
      acsg_curr ~ "ACS Global Current",
      acsg_retain ~ "Global Retained (%)",
      acsi_prev ~ "ACS IADL Before",
      acsi_curr ~ "ACS IADL Current",
      acsi_retain ~ "IADL Retained (%)",
      acsl_prev ~ "ACS Leisure Before",
      acsl_curr ~ "ACS Leisure Current",
      acsl_retain ~ "Leisure Retained (%)",
      acsf_prev ~ "ACS Fitness Before",
      acsf_curr ~ "ACS Fitness Current",
      acsf_retain ~ "Fitness Retained (%)",
      acss_prev ~ "ACS Social Before",
      acss_curr ~ "ACS Social Current",
      acss_retain ~ "Social Retained (%)"
    )
  ) 


print(mean_sd_acs)   

Characteristic	N = 451
ACS Global Before	70 (10)
ACS Global Current	53 (14)
Global Retained (%)	76 (17)
ACS IADL Before	21.51 (2.95)
ACS IADL Current	17.2 (4.7)
IADL Retained (%)	80 (18)
ACS Leisure Before	21.7 (5.1)
ACS Leisure Current	17.3 (5.4)
Leisure Retained (%)	80 (19)
ACS Fitness Before	13.1 (4.3)
ACS Fitness Current	8.3 (4.0)
Fitness Retained (%)	67 (33)
ACS Social Before	13.33 (1.52)
ACS Social Current	9.83 (2.74)
Social Retained (%)	74 (19)
1 Mean (SD)

Correlations

All variables

Correlation of all variables of interest with TBI QOL subscores. While too small to read in HTML print out, nice reference during analysis

library(ggplot2)
library(reshape2)
library(Hmisc)

all_variables <- c("age_current", "time_injury", "gender", "edu", "work_current", "income", "severity", "substance", "acsg_prev", "acsg_curr", "acsg_retain", "acsi_prev", "acsi_curr", "acsi_retain", "acsl_prev", "acsl_curr", "acsl_retain", "acsf_prev", "acsf_curr", "acsf_retain", "acss_prev", "acss_curr", "acss_retain","tbiqol_part_sra_tscore", "tbiqol_anxiety_tscore", "tbiqol_comm_tscore", "tbiqol_ue_tscore","tbiqol_depression_tscore", "tbiqol_fatigue_tscore", "tbiqol_genconcern_tscore", "tbiqol_grief_tscore", "tbiqol_mobility_tscore", "tbiqol_headache_tscore", "tbiqol_pain_tscore", "tbiqol_resilience_tscore", "tbiqol_satissra_tscore", "tbiqol_selfesteem_tscore", "tbiqol_stigma_tscore", "spstotal", "bfi_extraversion", "bfi_agreeable", "bfi_consciousness", "bfi_neuroticism", "bfi_openness", "frsbe_apathy", "frsbe_exec", "frsbe_disinhib", "frsbe_total")

all_variables <- intersect(all_variables, colnames(df))  # Ensure selected variables are in the dataframe

all_variables_df <- df[, all_variables]

cor_matrix <- rcorr(as.matrix(all_variables_df),type = "spearman")$r
cor_matrix[upper.tri(cor_matrix)] <- NA
p_matrix <- rcorr(as.matrix(all_variables_df),type = "spearman")$P
p_matrix[is.na(p_matrix)]=.0000001
p_matrix[upper.tri(p_matrix)] <- NA

# Melt the correlation matrix for ggplot
melted_cor <- melt(cor_matrix,na.rm = T)
melted_p = melt(p_matrix,na.rm = T)

melted_cor$p=melted_p$value
melted_cor$psig=""
melted_cor$psig[melted_cor$p<.05]="*"
melted_cor$psig[melted_cor$p<.01]="**"
melted_cor$psig[melted_cor$p<.001]="***"

# Create a heatmap using ggplot2
ggplot(melted_cor, aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(value,2)), vjust = 1, size = 2) +  # Adjust size here
  geom_text(aes(label = psig), vjust = .25, size = 6) +  # Adjust size here
  scale_fill_gradient2(low = "purple", mid = "white", high = "orange", 
                       midpoint = 0, limit = c(-1,1), space = "Lab",
                       name = "Correlation") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),  # Adjust size here
        axis.text.y = element_text(angle = 45, hjust = 1, size = 10),  # Adjust size here
        axis.title = element_text(size = 14),
        axis.ticks = element_line(size = 1)) +
  labs(caption = "<.05 = *, <.01 = **, <.001 = ***") +
  xlab("") +
  ylab("") +
  ggtitle("Correlation Matrix all variables")

This is pretty hard to read, so the following matrices break it down into smaller parts. First, we’ll rename variables for ease of reading

#Rename Variables in new dataset
df2 <- read_excel('cleaned_data_Final.xlsx')

# Rename multiple variables
df2 <- df2 %>%
  rename(Global = acsg_retain,
         IADL = acsi_retain,
         Leisure = acsl_retain,
         Fitness = acsf_retain,
         Social = acss_retain,
         Age = age_current,
         TimeSinceInjury = time_injury,
         Extraversion = bfi_extraversion,
         Agreeable = bfi_agreeable,
         Consciousness =bfi_consciousness,
         Neuroticism = bfi_neuroticism,
         Openness = bfi_openness,
         Apathy = frsbe_apathy,
         ExecFunc = frsbe_exec,
         Disinhibition = frsbe_disinhib,
         Total = frsbe_total,
         SocialSupport = spstotal,
         Communication = tbiqol_comm_tscore,
         ExecFuncQOL = tbiqol_execfunc_tscore, 
         GeneralCognition = tbiqol_genconcern_tscore,
         UpperExtremity = tbiqol_ue_tscore, 
         Fatigue = tbiqol_fatigue_tscore, 
         Mobility = tbiqol_mobility_tscore, 
         Headache = tbiqol_headache_tscore,
         Pain = tbiqol_pain_tscore,
         Anger = tbiqol_anger_tscore,
         PositiveAffect = tbiqol_posaffect_tscore,
         Age = age_current,
         Education = edu,
         Work = work_current,
         SubstanceUse = substance,
         Anxiety = tbiqol_anxiety_tscore, 
         Depression = tbiqol_depression_tscore, 
         Grief = tbiqol_grief_tscore, 
         TraitResilience = tbiqol_resilience_tscore,  
         SelfEsteem = tbiqol_selfesteem_tscore, 
         Stigma = tbiqol_stigma_tscore,
         TimeSinceInjury = time_injury,
         MaritalStatus = marital_status,
         SocialSupport = spstotal,
         HouseholdSize = house_size,
         PhysicalHealth = phys_health_index,
         EmotionalHealth = emo_health_index)

Included Variables

Matrix with heat map for all included variables in dissertation. Figure 4 in dissertation

library(ggplot2)
library(reshape2)
library(Hmisc)

# Define QOL_PPF_variables
QOL_PPF_variables <- c("Global", "IADL", "Leisure", "Fitness", "Social", "PhysicalHealth", "GeneralCognition", "EmotionalHealth", "Extraversion", "Agreeable", "Consciousness", "Neuroticism", "Openness", "income", "MaritalStatus", "SocialSupport", "Apathy", "Disinhibition", "ExecFunc", "Total")

# Ensure selected variables are in the dataframe
QOL_PPF_variables <- intersect(QOL_PPF_variables, colnames(df2))

# Extract relevant data
QOL_PPF_df2 <- df2[, QOL_PPF_variables]

# Calculate correlation matrix
cor_matrix <- rcorr(as.matrix(QOL_PPF_df2), type = "spearman")$r
cor_matrix[upper.tri(cor_matrix)] <- NA
p_matrix <- rcorr(as.matrix(QOL_PPF_df2), type = "spearman")$P
p_matrix[is.na(p_matrix)] <- .0000001
p_matrix[upper.tri(p_matrix)] <- NA

# Melt the correlation matrix for ggplot
melted_cor <- melt(cor_matrix, na.rm = TRUE)
melted_p <- melt(p_matrix, na.rm = TRUE)

melted_cor$p <- melted_p$value
melted_cor$psig <- ""
melted_cor$psig[melted_cor$p < .05] <- "*"
melted_cor$psig[melted_cor$p < .01] <- "**"
melted_cor$psig[melted_cor$p < .001] <- "***"

# Create a heatmap using ggplot2
ggplot(melted_cor, aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(value,2)), vjust = 1, size = 6) +  # Adjust size here
  geom_text(aes(label = psig), vjust = .25, size = 6) +  # Adjust size here
  scale_fill_gradient2(low = "purple", mid = "white", high = "orange", 
                       midpoint = 0, limit = c(-1,1), space = "Lab",
                       name = "Correlation") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 20),  # Adjust size here
        axis.text.y = element_text(angle = 45, hjust = 1, size = 20),  # Adjust size here
        axis.title = element_text(size = 14),
        axis.ticks = element_line(size = 1)) +
  labs(caption = "<.05 = *, <.01 = **, <.001 = ***") +
  xlab("") +
  ylab("") +
  ggtitle("Correlation Matrix Personal Protective Factors with ACS Variables") +
  guides(fill = guide_colorbar(title.position = "right"))

Below is the breakdown of all TBIQOL sub scores with the ACS3. While not included in this study, helpful for discussion and future publications.

# Define QOL_variables
QOL_variables <- c("Global", "IADL", "Leisure", "Fitness", "Social", "Anxiety", "Depression", "Grief", "Resilience", "SelfEsteem", "Stigma", "Communication", "GeneralCognition", "ExecFuncQOL", "UpperExtremity", "Fatigue", "Mobility", "Headache", "Pain")

# Ensure selected variables are in the dataframe
QOL_variables <- intersect(QOL_variables, colnames(df2))

# Extract relevant data
QOL_df2 <- df2[, QOL_variables]

# Calculate correlation matrix
cor_matrix <- rcorr(as.matrix(QOL_df2), type = "spearman")$r
cor_matrix[upper.tri(cor_matrix)] <- NA
p_matrix <- rcorr(as.matrix(QOL_df2), type = "spearman")$P
p_matrix[is.na(p_matrix)] <- .0000001
p_matrix[upper.tri(p_matrix)] <- NA

# Melt the correlation matrix for ggplot
melted_cor <- melt(cor_matrix, na.rm = TRUE)
melted_p <- melt(p_matrix, na.rm = TRUE)

melted_cor$p <- melted_p$value
melted_cor$psig <- ""
melted_cor$psig[melted_cor$p < .05] <- "*"
melted_cor$psig[melted_cor$p < .01] <- "**"
melted_cor$psig[melted_cor$p < .001] <- "***"

# Create a heatmap using ggplot2
ggplot(melted_cor, aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(value,2)), vjust = 1, size = 6) +  # Adjust size here
  geom_text(aes(label = psig), vjust = .25, size = 6) +  # Adjust size here
  scale_fill_gradient2(low = "purple", mid = "white", high = "orange", 
                       midpoint = 0, limit = c(-1,1), space = "Lab",
                       name = "Correlation") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 20),  # Adjust size here
        axis.text.y = element_text(angle = 45, hjust = 1, size = 20),  # Adjust size here
        axis.title = element_text(size = 14),
        axis.ticks = element_line(size = 1)) +
  labs(caption = "<.05 = *, <.01 = **, <.001 = ***") +
  xlab("") +
  ylab("") +
  ggtitle("Correlation Matrix Personal Protective Factors with ACS Variables")

FrSBe and TBI-QOL

TBI QOL subscales with the FrSBe. I think this could be really interested as a future paper given how some subscales of the TBI QOL overlap with FrSBe (Exec functioning, anger, dyscontrol)

QOL_FRSBE_variables <- c("frsbe_apathy", "frsbe_disinhib", "frsbe_exec", "frsbe_total", "tbiqol_part_sra_tscore", "tbiqol_anxiety_tscore", "tbiqol_comm_tscore","tbiqol_ue_tscore","tbiqol_depression_tscore", "tbiqol_dyscontrol_tscoree", "tbiqol_execfunc_tscore", "tbiqol_fatigue_tscore", "tbiqol_genconcern_tscore", "tbiqol_grief_tscore", "tbiqol_mobility_tscore", "tbiqol_headache_tscore", "tbiqol_pain_tscore", "tbiqol_resilience_tscore", "tbiqol_satissra_tscore", "tbiqol_selfesteem_tscore", "tbiqol_stigma_tscore")

QOL_FRSBE_variables <- intersect(QOL_FRSBE_variables, colnames(df))  # Ensure selected variables are in the dataframe

QOL_FRSBE_df <- df[, QOL_FRSBE_variables]

cor_matrix <- rcorr(as.matrix(QOL_FRSBE_df ),type = "spearman")$r
cor_matrix[upper.tri(cor_matrix)] <- NA
p_matrix <- rcorr(as.matrix(QOL_FRSBE_df ),type = "spearman")$P
p_matrix[is.na(p_matrix)]=.0000001
p_matrix[upper.tri(p_matrix)] <- NA

# Melt the correlation matrix for ggplot
melted_cor <- melt(cor_matrix,na.rm = T)
melted_p = melt(p_matrix,na.rm = T)

melted_cor$p=melted_p$value
melted_cor$psig=""
melted_cor$psig[melted_cor$p<.05]="*"
melted_cor$psig[melted_cor$p<.01]="**"
melted_cor$psig[melted_cor$p<.001]="***"

# Create a heatmap using ggplot2
ggplot(melted_cor, aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(value,2)), vjust = 1) +
  geom_text(aes(label = psig), vjust = .25,size=5) +
  scale_fill_gradient2(low = "purple", mid = "white", high = "orange", 
                       midpoint = 0, limit = c(-1,1), space = "Lab",
                       name="Correlation") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_text(angle = 45, hjust = 1))+
  labs(caption = "<.05 = *, <.01 = **, <.001 = ***")+
  xlab("")+
  ylab("")+
  ggtitle("Correlation Matrix TBI-QOL Subscale and FrSBe Variables")

RQ1: Regression Analysis

Research Question 1 1. What is the relationship between protective factors and self-regulation with resiliency-related outcomes such as re-engagement in meaningful activities? a. To what extent do protective factors and self-regulation predict resiliency-related outcomes in the TBI population? Hypothesis: Higher self-regulation will be associated with better resiliency-related outcomes b. To what extent does self-regulation mediate or moderate the influence of protective factors on resiliency-related outcomes after TBI? Hypothesis: Self-regulation will impact the relationship between protective factors and resiliency-related outcomes

RQ1a

First, we’ll look at the hierarchical linear model as outlined in Chapter 3. Then, to dive deeper, a “post hoc” analysis of each subscale of the ACS and use AIC to determine model of best fit.

ACS Global

In this section, we’ll do the original hierarchical model with protective and environmental protective factors in the first step and then total self-regulation score added for the second. *note that the cognitive composite score is not included as it includes exec functioning, which in this paper is seen as a self-regulatory process. therefore, gen concerns (memory and concentration) is used as cognitive protective factor

step1 <- lm(acsg_retain~age_current, data=df)
summary(step1)

## 
## Call:
## lm(formula = acsg_retain ~ age_current, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.534 -11.401  -0.019   7.551  36.793 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  92.4076     8.1805  11.296 1.88e-14 ***
## age_current  -0.3515     0.1656  -2.122   0.0396 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.98 on 43 degrees of freedom
## Multiple R-squared:  0.0948, Adjusted R-squared:  0.07375 
## F-statistic: 4.504 on 1 and 43 DF,  p-value: 0.03962

step2<- lm(acsg_retain~age_current+phys_health_index+tbiqol_genconcern_tscore+emo_health_index+ spstotal, data=df)
summary(step2)

## 
## Call:
## lm(formula = acsg_retain ~ age_current + phys_health_index + 
##     tbiqol_genconcern_tscore + emo_health_index + spstotal, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.222 -10.023  -0.194   9.987  31.974 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)              36.64130   25.97605   1.411   0.1663  
## age_current              -0.08347    0.17014  -0.491   0.6264  
## phys_health_index         0.04845    0.22846   0.212   0.8332  
## tbiqol_genconcern_tscore  0.94332    0.38078   2.477   0.0177 *
## emo_health_index         -0.06226    0.22244  -0.280   0.7810  
## spstotal                  0.12949    0.24409   0.531   0.5988  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.52 on 39 degrees of freedom
## Multiple R-squared:  0.3224, Adjusted R-squared:  0.2355 
## F-statistic: 3.711 on 5 and 39 DF,  p-value: 0.00764

#Nested Model Comparison
anova(step1, step2)

## Analysis of Variance Table
## 
## Model 1: acsg_retain ~ age_current
## Model 2: acsg_retain ~ age_current + phys_health_index + tbiqol_genconcern_tscore + 
##     emo_health_index + spstotal
##   Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
## 1     43 10984.7                              
## 2     39  8222.7  4    2762.1 3.2751 0.02082 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#change in R-squared
summary(step2)$r.squared - summary(step1)$r.squared

## [1] 0.2276085

step3<- lm(acsg_retain~age_current+phys_health_index+tbiqol_genconcern_tscore+emo_health_index+ spstotal+frsbe_total, data=df)
summary(step3)

## 
## Call:
## lm(formula = acsg_retain ~ age_current + phys_health_index + 
##     tbiqol_genconcern_tscore + emo_health_index + spstotal + 
##     frsbe_total, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.897  -9.526   0.100   8.754  32.154 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)              53.31446   38.64125   1.380   0.1757  
## age_current              -0.08802    0.17176  -0.512   0.6113  
## phys_health_index         0.05187    0.23048   0.225   0.8231  
## tbiqol_genconcern_tscore  0.87069    0.40347   2.158   0.0373 *
## emo_health_index         -0.08405    0.22738  -0.370   0.7137  
## spstotal                  0.09315    0.25383   0.367   0.7157  
## frsbe_total              -0.08471    0.14432  -0.587   0.5607  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.64 on 38 degrees of freedom
## Multiple R-squared:  0.3285, Adjusted R-squared:  0.2225 
## F-statistic: 3.098 on 6 and 38 DF,  p-value: 0.01436

#Nested Model Comparison
anova(step2, step3)

## Analysis of Variance Table
## 
## Model 1: acsg_retain ~ age_current + phys_health_index + tbiqol_genconcern_tscore + 
##     emo_health_index + spstotal
## Model 2: acsg_retain ~ age_current + phys_health_index + tbiqol_genconcern_tscore + 
##     emo_health_index + spstotal + frsbe_total
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     39 8222.7                           
## 2     38 8148.8  1    73.884 0.3445 0.5607

#change in R-squared
summary(step3)$r.squared - summary(step2)$r.squared

## [1] 0.006088374

Assumptions tested

Test for multicollinearity

library(car)
vif(step3)

##              age_current        phys_health_index tbiqol_genconcern_tscore 
##                 1.280847                 2.110368                 2.477995 
##         emo_health_index                 spstotal              frsbe_total 
##                 2.090700                 1.712354                 1.670229

Post Hoc AIC Models

As we have a smaller n and need to be parsimonious with the variables we use in the regression model, we’ll look at several models based on correlations (higher correlations added to the models) and then calculate the AIC. The lower AIC, the better the fit and that model will be used

Note that all assumptions were tested for each model and were met

ACS3 Social

library(AICcmodavg)
#composite scores
model1 <- lm(acss_retain~phys_health_index+cog_health_index+emo_health_index+spstotal, data=df)

model2 <- lm(acss_retain~phys_health_index+cog_health_index+emo_health_index+spstotal+frsbe_total, data=df)

model3 <- lm(acss_retain~phys_health_index+cog_health_index+emo_health_index+spstotal+frsbe_apathy, data=df)

#individual
model4 <- lm(acss_retain~tbiqol_fatigue_tscore+tbiqol_genconcern_tscore+tbiqol_depression_tscore+tbiqol_anxiety_tscore+spstotal+spstotal, data=df)
model5 <- lm(acss_retain~tbiqol_fatigue_tscore+tbiqol_genconcern_tscore+tbiqol_depression_tscore+tbiqol_anxiety_tscore+spstotal+frsbe_total, data=df)
model6 <- lm(acss_retain~tbiqol_fatigue_tscore+tbiqol_genconcern_tscore+tbiqol_depression_tscore+tbiqol_anxiety_tscore+spstotal+frsbe_apathy, data=df)
model7 <- lm(acss_retain~tbiqol_fatigue_tscore+tbiqol_genconcern_tscore+tbiqol_depression_tscore+tbiqol_anxiety_tscore+frsbe_apathy, data=df)
model8 <- lm(acss_retain~tbiqol_fatigue_tscore+tbiqol_genconcern_tscore+tbiqol_depression_tscore+tbiqol_anxiety_tscore+frsbe_total, data=df) 
model9 <- lm(acss_retain~tbiqol_fatigue_tscore+tbiqol_comm_tscore+tbiqol_depression_tscore+tbiqol_anxiety_tscore+frsbe_total, data=df)
model10 <- lm(acss_retain~tbiqol_fatigue_tscore+tbiqol_comm_tscore+tbiqol_depression_tscore+tbiqol_anxiety_tscore+frsbe_total, data=df)
model11 <- lm(acss_retain~tbiqol_fatigue_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+ tbiqol_depression_tscore+tbiqol_anxiety_tscore+frsbe_apathy, data=df)

models <- list (model1, model2, model3, model4, model5, model6, model7, model8, model9, model10, model11)

mod.names <- c('model1', 'model2', 'model3', 'model4', 'model5','model6', 'model7', 'model8', 'model9', 'model10', 'model11')

aictab(cand.set = models, modnames = mod.names)

## 
## Model selection based on AICc:
## 
##         K   AICc Delta_AICc AICcWt Cum.Wt      LL
## model7  7 381.10       0.00   0.44   0.44 -182.04
## model11 8 382.14       1.04   0.26   0.70 -181.07
## model6  8 384.02       2.92   0.10   0.80 -182.01
## model9  7 385.44       4.34   0.05   0.85 -184.20
## model10 7 385.44       4.34   0.05   0.90 -184.20
## model3  7 386.05       4.95   0.04   0.94 -184.51
## model8  7 387.20       6.10   0.02   0.96 -185.08
## model4  7 387.60       6.50   0.02   0.97 -185.29
## model1  6 387.71       6.61   0.02   0.99 -186.75
## model5  8 390.06       8.96   0.00   1.00 -185.03
## model2  7 390.10       9.00   0.00   1.00 -186.54

summary(model7)

## 
## Call:
## lm(formula = acss_retain ~ tbiqol_fatigue_tscore + tbiqol_genconcern_tscore + 
##     tbiqol_depression_tscore + tbiqol_anxiety_tscore + frsbe_apathy, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.469  -9.018  -0.100  12.942  24.694 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)              100.94678   30.60729   3.298  0.00208 **
## tbiqol_fatigue_tscore      0.06575    0.36559   0.180  0.85820   
## tbiqol_genconcern_tscore   0.68936    0.33991   2.028  0.04942 * 
## tbiqol_depression_tscore   0.37274    0.34962   1.066  0.29292   
## tbiqol_anxiety_tscore     -0.87131    0.34569  -2.520  0.01592 * 
## frsbe_apathy              -0.83974    0.33446  -2.511  0.01630 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.85 on 39 degrees of freedom
## Multiple R-squared:  0.4592, Adjusted R-squared:  0.3898 
## F-statistic: 6.622 on 5 and 39 DF,  p-value: 0.0001501

ACS3 IADL

library(AICcmodavg)
#composite scores
model1 <- lm(acsi_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+tbiqol_grief_tscore+spstotal+frsbe_total, data=df)
model2 <- lm(acsi_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+tbiqol_grief_tscore+spstotal+frsbe_apathy, data=df)
model3<- lm(acsi_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+tbiqol_grief_tscore+spstotal+frsbe_exec, data=df)
model4 <- lm(acsi_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+spstotal+tbiqol_ue_tscore+frsbe_apathy, data=df)
model5 <- lm(acsi_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+tbiqol_ue_tscore+spstotal+frsbe_apathy, data=df)
model6<- lm(acsi_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+tbiqol_ue_tscore+frsbe_apathy, data=df)

models <- list (model1, model2, model3, model4, model5, model6)

mod.names <- c('model1', 'model2', 'model3', 'model4', 'model5','model6')

aictab(cand.set = models, modnames = mod.names)

## 
## Model selection based on AICc:
## 
##        K   AICc Delta_AICc AICcWt Cum.Wt      LL
## model6 7 385.38       0.00   0.42   0.42 -184.18
## model4 7 386.02       0.64   0.30   0.72 -184.50
## model5 8 388.35       2.97   0.09   0.82 -184.17
## model3 8 388.74       3.36   0.08   0.90 -184.37
## model2 8 389.19       3.81   0.06   0.96 -184.60
## model1 8 390.01       4.63   0.04   1.00 -185.01

summary(model6)

## 
## Call:
## lm(formula = acsi_retain ~ tbiqol_mobility_tscore + tbiqol_genconcern_tscore + 
##     tbiqol_comm_tscore + tbiqol_ue_tscore + frsbe_apathy, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -35.730 -10.279  -0.783  11.626  25.604 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)               41.1421    22.2568   1.849   0.0721 .
## tbiqol_mobility_tscore     0.3760     0.4477   0.840   0.4061  
## tbiqol_genconcern_tscore   0.7421     0.3466   2.141   0.0386 *
## tbiqol_comm_tscore        -0.2813     0.3751  -0.750   0.4579  
## tbiqol_ue_tscore           0.4750     0.4531   1.048   0.3009  
## frsbe_apathy              -0.3969     0.3296  -1.204   0.2358  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.57 on 39 degrees of freedom
## Multiple R-squared:  0.3593, Adjusted R-squared:  0.2772 
## F-statistic: 4.374 on 5 and 39 DF,  p-value: 0.002952

ACS3 Leisure

model1 <- lm (acsl_retain~tbiqol_mobility_tscore+tbiqol_ue_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+frsbe_total, data=df)
model2 <- lm (acsl_retain~tbiqol_mobility_tscore+tbiqol_ue_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+tbiqol_grief_tscore+frsbe_total, data=df)
model3 <-lm (acsl_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+frsbe_total, data=df)
model4 <-lm (acsl_retain~tbiqol_mobility_tscore+tbiqol_grief_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+frsbe_total, data=df)
model5 <-lm (acsl_retain~tbiqol_ue_tscore+tbiqol_grief_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+frsbe_total, data=df)
model6 <-lm (acsl_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+frsbe_apathy, data=df)
model7 <-lm (acsl_retain~tbiqol_mobility_tscore+tbiqol_genconcern_tscore+tbiqol_comm_tscore+frsbe_exec, data=df)

models <- list (model1, model2, model3, model4, model5, model6, model7)

mod.names <- c('model1', 'model2', 'model3', 'model4', 'model5', 'model6', 'model7')

aictab(cand.set = models, modnames = mod.names)

## 
## Model selection based on AICc:
## 
##        K   AICc Delta_AICc AICcWt Cum.Wt      LL
## model6 6 389.41       0.00   0.31   0.31 -187.60
## model7 6 389.72       0.30   0.27   0.58 -187.75
## model3 6 389.98       0.56   0.24   0.82 -187.88
## model1 7 392.16       2.75   0.08   0.89 -187.57
## model4 7 392.34       2.92   0.07   0.97 -187.65
## model2 8 394.90       5.48   0.02   0.99 -187.45
## model5 7 395.80       6.39   0.01   1.00 -189.39

summary(model6)

## 
## Call:
## lm(formula = acsl_retain ~ tbiqol_mobility_tscore + tbiqol_genconcern_tscore + 
##     tbiqol_comm_tscore + frsbe_apathy, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.380 -10.329  -1.113  11.003  36.264 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)               45.7099    23.7147   1.927   0.0610 .
## tbiqol_mobility_tscore     0.6242     0.3538   1.764   0.0853 .
## tbiqol_genconcern_tscore   0.8329     0.3692   2.256   0.0296 *
## tbiqol_comm_tscore        -0.3214     0.3869  -0.831   0.4111  
## frsbe_apathy              -0.2811     0.3454  -0.814   0.4206  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.59 on 40 degrees of freedom
## Multiple R-squared:  0.2816, Adjusted R-squared:  0.2098 
## F-statistic:  3.92 on 4 and 40 DF,  p-value: 0.008891

ACS3 Fitness

model1 <- lm (acsf_retain~tbiqol_mobility_tscore+tbiqol_pain_tscore+tbiqol_genconcern_tscore+tbiqol_anxiety_tscore+tbiqol_stigma_tscore+frsbe_total, data=df)
model2 <- lm (acsf_retain~phys_health_index+tbiqol_genconcern_tscore+tbiqol_anxiety_tscore+tbiqol_stigma_tscore, data=df)
model3 <-lm (acsf_retain~tbiqol_mobility_tscore+tbiqol_pain_tscore+tbiqol_genconcern_tscore+tbiqol_anxiety_tscore+tbiqol_grief_tscore, data=df)
model4 <-lm (acsf_retain~phys_health_index+tbiqol_genconcern_tscore+tbiqol_anxiety_tscore+tbiqol_grief_tscore, data=df)
model5 <- lm (acsf_retain~tbiqol_mobility_tscore+tbiqol_pain_tscore+tbiqol_genconcern_tscore+tbiqol_anxiety_tscore+tbiqol_grief_tscore+tbiqol_stigma_tscore, data=df)
model6 <- lm (acsf_retain~tbiqol_mobility_tscore+tbiqol_pain_tscore+tbiqol_genconcern_tscore+tbiqol_anxiety_tscore+tbiqol_stigma_tscore+frsbe_total+income, data=df)
model7 <- lm (acsf_retain~phys_health_index+tbiqol_genconcern_tscore+tbiqol_anxiety_tscore+tbiqol_stigma_tscore+frsbe_total+income, data=df)
model8 <- lm (acsf_retain~phys_health_index+tbiqol_genconcern_tscore+tbiqol_anxiety_tscore+tbiqol_stigma_tscore+income, data=df)

models <- list (model1, model2, model3, model4, model5, model6, model7, model8)

mod.names <- c('model1', 'model2', 'model3', 'model4', 'model5', 'model6', 'model7', 'model8')

aictab(cand.set = models, modnames = mod.names)

## 
## Model selection based on AICc:
## 
##        K   AICc Delta_AICc AICcWt Cum.Wt      LL
## model8 7 440.01       0.00   0.46   0.46 -211.45
## model2 6 440.66       0.65   0.33   0.79 -213.19
## model7 8 442.90       2.89   0.11   0.90 -211.39
## model1 8 445.26       5.25   0.03   0.94 -212.57
## model5 8 445.51       5.50   0.03   0.97 -212.70
## model6 9 445.57       5.56   0.03   0.99 -211.14
## model4 6 449.75       9.74   0.00   1.00 -217.77
## model3 7 451.25      11.24   0.00   1.00 -217.11

summary(model8)

## 
## Call:
## lm(formula = acsf_retain ~ phys_health_index + tbiqol_genconcern_tscore + 
##     tbiqol_anxiety_tscore + tbiqol_stigma_tscore + income, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -50.483 -15.525  -6.805   5.882 123.848 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)               49.0461    68.8141   0.713   0.4804  
## phys_health_index          0.2584     0.4671   0.553   0.5833  
## tbiqol_genconcern_tscore   0.4865     0.7595   0.641   0.5257  
## tbiqol_anxiety_tscore     -0.1576     0.7014  -0.225   0.8234  
## tbiqol_stigma_tscore      -0.7559     0.8658  -0.873   0.3881  
## income                    12.1671     6.8708   1.771   0.0846 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.82 on 38 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.216,  Adjusted R-squared:  0.1128 
## F-statistic: 2.094 on 5 and 38 DF,  p-value: 0.08746

RQ1b

Moderation

For moderation, looking at personal protective factors and environmental protective factors from original model

mod_acsg_phys <- lm(acsg_retain~phys_health_index*frsbe_total, data=df)
summary(mod_acsg_phys)

## 
## Call:
## lm(formula = acsg_retain ~ phys_health_index * frsbe_total, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.220  -8.719  -1.826  10.001  34.734 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)                    47.828275 105.856796   0.452    0.654
## phys_health_index               0.554929   1.052087   0.527    0.601
## frsbe_total                    -0.031261   0.948658  -0.033    0.974
## phys_health_index:frsbe_total  -0.002117   0.009561  -0.221    0.826
## 
## Residual standard error: 15.28 on 41 degrees of freedom
## Multiple R-squared:  0.2115, Adjusted R-squared:  0.1538 
## F-statistic: 3.666 on 3 and 41 DF,  p-value: 0.01982

mod_acsg_genconcern <- lm(acsg_retain~tbiqol_genconcern_tscore*frsbe_total, data=df)
summary(mod_acsg_genconcern)

## 
## Call:
## lm(formula = acsg_retain ~ tbiqol_genconcern_tscore * frsbe_total, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.201  -8.855   0.300   8.808  33.111 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)
## (Intercept)                          81.233098  55.408551   1.466    0.150
## tbiqol_genconcern_tscore              0.089380   1.453353   0.061    0.951
## frsbe_total                          -0.392674   0.511644  -0.767    0.447
## tbiqol_genconcern_tscore:frsbe_total  0.008742   0.014161   0.617    0.540
## 
## Residual standard error: 14.12 on 41 degrees of freedom
## Multiple R-squared:  0.3265, Adjusted R-squared:  0.2772 
## F-statistic: 6.625 on 3 and 41 DF,  p-value: 0.0009409

mod_acsg_emo <- lm(acsg_retain~emo_health_index*frsbe_total, data=df)
summary(mod_acsg_emo)

## 
## Call:
## lm(formula = acsg_retain ~ emo_health_index * frsbe_total, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.029  -7.277   0.207   7.957  30.884 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)
## (Intercept)                  -46.603657 116.746539  -0.399    0.692
## emo_health_index               1.416602   1.080829   1.311    0.197
## frsbe_total                    0.903860   1.018454   0.887    0.380
## emo_health_index:frsbe_total  -0.010967   0.009591  -1.143    0.259
## 
## Residual standard error: 15.46 on 41 degrees of freedom
## Multiple R-squared:  0.1929, Adjusted R-squared:  0.1339 
## F-statistic: 3.267 on 3 and 41 DF,  p-value: 0.03075

mod_acsg_sps <- lm(acsg_retain~spstotal*frsbe_total, data=df)
summary(mod_acsg_sps)

## 
## Call:
## lm(formula = acsg_retain ~ spstotal * frsbe_total, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.883  -9.235   0.023   8.489  30.764 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)           49.500184 115.410475   0.429    0.670
## spstotal               0.670089   1.422480   0.471    0.640
## frsbe_total            0.067652   1.044274   0.065    0.949
## spstotal:frsbe_total  -0.004072   0.013179  -0.309    0.759
## 
## Residual standard error: 15.7 on 41 degrees of freedom
## Multiple R-squared:  0.1672, Adjusted R-squared:  0.1063 
## F-statistic: 2.745 on 3 and 41 DF,  p-value: 0.05523

There was no moderating effect of apathy on any of the predictors.

Mediation

Here we look at mediation effect of the total FrSBe scores on the personal and environmental factors used in the post hoc model (mobility, general cog functioning, anxiety, depression, and social support)

How to read: ACME = indirect effect. ADE = direct effect. ACME + ADE = total effect.

Physical Health

# Initial Model
model1 <- lm(acsg_retain ~ phys_health_index, df) # Y ~ X, DV predicted by IV - no mediation considered - total effect
summary(model1)

## 
## Call:
## lm(formula = acsg_retain ~ phys_health_index, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.889  -9.129  -2.889  10.208  34.025 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        33.5455    16.1736   2.074   0.0441 *
## phys_health_index   0.4457     0.1688   2.640   0.0115 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.58 on 43 degrees of freedom
## Multiple R-squared:  0.1395, Adjusted R-squared:  0.1195 
## F-statistic: 6.969 on 1 and 43 DF,  p-value: 0.01151

# Mediation paths
medmodel1 <- lm(frsbe_total ~ phys_health_index, df) # M ~ X, mediator predicted by X
outputmodel1 <- lm(acsg_retain ~ phys_health_index + frsbe_total, df) # Y ~ X + M, DV predicted by mediator, adjusting for IV

summary(medmodel1)

## 
## Call:
## lm(formula = frsbe_total ~ phys_health_index, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.474 -13.988  -1.927  11.537  48.542 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       154.7515    19.4114   7.972 5.19e-10 ***
## phys_health_index  -0.5028     0.2026  -2.481   0.0171 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.7 on 43 degrees of freedom
## Multiple R-squared:  0.1252, Adjusted R-squared:  0.1049 
## F-statistic: 6.156 on 1 and 43 DF,  p-value: 0.01709

summary(outputmodel1)

## 
## Call:
## lm(formula = acsg_retain ~ phys_health_index + frsbe_total, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.407  -8.551  -1.666   9.717  35.212 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)        70.6031    24.6747   2.861  0.00655 **
## phys_health_index   0.3253     0.1750   1.860  0.06997 . 
## frsbe_total        -0.2395     0.1231  -1.945  0.05854 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.1 on 42 degrees of freedom
## Multiple R-squared:  0.2106, Adjusted R-squared:  0.173 
## F-statistic: 5.601 on 2 and 42 DF,  p-value: 0.006979

# Mediation test
mediation <- mediate(medmodel1, # Mediator model
                    outputmodel1, # Outcome model
                    boot = T, # Ask for bootstrapped confidence intervals
                    treat="phys_health_index", # Name of the x variable
                    mediator="frsbe_total" # Name of the m variable
                    )
# if you don't want bootstrap, just delete 'sims' line and set boot = F

summary(mediation)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value  
## ACME            0.12040      0.01639         0.30   0.016 *
## ADE             0.32533     -0.00944         0.66   0.062 .
## Total Effect    0.44572      0.10061         0.80   0.016 *
## Prop. Mediated  0.27011      0.03183         0.95   0.032 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 45 
## 
## 
## Simulations: 1000

plot(mediation)

There is a significant indirect effect and an insignificant direct effect, indicating total mediation

library(diagram)

patha <- paste(round(as.numeric(coef(medmodel1)[2]),2))
patha

## [1] "-0.5"

pathb <- paste(round(as.numeric(coef(outputmodel1)[3]),2))
pathb

## [1] "-0.24"

# indirect
direct <- paste(round(mediation$z.avg,2)) # direct
total <- paste(round((mediation$d.avg + mediation$z.avg),2)) # total

pathc <- paste(direct,"(",total,")")

data <- c(0, patha, 0,
          0, 0, 0, 
          pathb, pathc, 0) # order: top = effect of IV on mediator; bottom left = effect of mediator on DV; bottom middle = direct effect (total effect)
M<- matrix (nrow=3, ncol=3, byrow = TRUE, data=data)
plot<- plotmat (M, pos=c(1,2), 
                name= c("Self-regulation: Total","Physical Health", "ACS Global"), # order: Mediator, IV, DV 
                box.type = "rect", box.size = 0.12, box.prop=0.5,  curve=0)

#### Cognitive Health

# MEDIATION
# Initial Model
model2 <- lm(acsg_retain ~ tbiqol_genconcern_tscore, df) # Y ~ X, DV predicted by IV - no mediation considered - total effect
summary(model2)

## 
## Call:
## lm(formula = acsg_retain ~ tbiqol_genconcern_tscore, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.388 -10.545   0.593   8.741  32.307 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               36.4502     9.1336   3.991 0.000252 ***
## tbiqol_genconcern_tscore   1.0783     0.2437   4.424 6.51e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.93 on 43 degrees of freedom
## Multiple R-squared:  0.3128, Adjusted R-squared:  0.2968 
## F-statistic: 19.57 on 1 and 43 DF,  p-value: 6.506e-05

# Mediation paths
medmodel2 <- lm(frsbe_total ~ tbiqol_genconcern_tscore, df) # M ~ X, mediator predicted by X
outputmodel2 <- lm(acsg_retain ~ tbiqol_genconcern_tscore + frsbe_total, df) # Y ~ X + M, DV predicted by mediator, adjusting for IV

summary(medmodel2)

## 
## Call:
## lm(formula = frsbe_total ~ tbiqol_genconcern_tscore, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.277 -11.050   0.906  10.408  41.270 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              153.4595    10.9219   14.05  < 2e-16 ***
## tbiqol_genconcern_tscore  -1.2707     0.2915   -4.36 7.98e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.65 on 43 degrees of freedom
## Multiple R-squared:  0.3065, Adjusted R-squared:  0.2904 
## F-statistic: 19.01 on 1 and 43 DF,  p-value: 7.976e-05

summary(outputmodel2)

## 
## Call:
## lm(formula = acsg_retain ~ tbiqol_genconcern_tscore + frsbe_total, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.675  -9.565   0.826   7.693  32.876 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)   
## (Intercept)              49.81073   21.73399   2.292  0.02699 * 
## tbiqol_genconcern_tscore  0.96771    0.29456   3.285  0.00206 **
## frsbe_total              -0.08706    0.12834  -0.678  0.50125   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.01 on 42 degrees of freedom
## Multiple R-squared:  0.3202, Adjusted R-squared:  0.2879 
## F-statistic: 9.893 on 2 and 42 DF,  p-value: 0.0003016

# Mediation test
mediation <- mediate(medmodel2, # Mediator model
                    outputmodel2, # Outcome model
                    boot = T, # Ask for bootstrapped confidence intervals
                    treat="tbiqol_genconcern_tscore", # Name of the x variable
                    mediator="frsbe_total" # Name of the m variable
                    )
# if you don't want bootstrap, just delete 'sims' line and set boot = F

summary(mediation)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME              0.111       -0.107         0.37    0.27    
## ADE               0.968        0.516         1.48  <2e-16 ***
## Total Effect      1.078        0.644         1.57  <2e-16 ***
## Prop. Mediated    0.103       -0.115         0.33    0.27    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 45 
## 
## 
## Simulations: 1000

plot(mediation)

There is no significant indirect effect. No mediation

Emotional Health

# Initial Model
model1 <- lm(acsg_retain ~ emo_health_index, df) # Y ~ X, DV predicted by IV - no mediation considered - total effect
summary(model1)

## 
## Call:
## lm(formula = acsg_retain ~ emo_health_index, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -35.283  -9.291  -0.427   8.717  31.818 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       38.5048    17.2256   2.235   0.0306 *
## emo_health_index   0.3741     0.1711   2.186   0.0343 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.94 on 43 degrees of freedom
## Multiple R-squared:    0.1,  Adjusted R-squared:  0.07909 
## F-statistic: 4.779 on 1 and 43 DF,  p-value: 0.03431

# Mediation paths
medmodel1 <- lm(frsbe_total ~ emo_health_index, df) # M ~ X, mediator predicted by X
outputmodel1 <- lm(acsg_retain ~ emo_health_index + frsbe_total, df) # Y ~ X + M, DV predicted by mediator, adjusting for IV

summary(medmodel1)

## 
## Call:
## lm(formula = frsbe_total ~ emo_health_index, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.178 -12.703  -2.888   7.872  48.247 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      176.3722    18.7984   9.382 5.76e-12 ***
## emo_health_index  -0.6950     0.1868  -3.721 0.000571 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.39 on 43 degrees of freedom
## Multiple R-squared:  0.2436, Adjusted R-squared:  0.226 
## F-statistic: 13.85 on 1 and 43 DF,  p-value: 0.0005709

summary(outputmodel1)

## 
## Call:
## lm(formula = acsg_retain ~ emo_health_index + frsbe_total, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.185  -7.578  -0.481   8.659  32.473 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)       82.6619    29.2672   2.824  0.00722 **
## emo_health_index   0.2001     0.1915   1.045  0.30208   
## frsbe_total       -0.2504     0.1360  -1.841  0.07273 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.51 on 42 degrees of freedom
## Multiple R-squared:  0.1672, Adjusted R-squared:  0.1275 
## F-statistic: 4.216 on 2 and 42 DF,  p-value: 0.02145

# Mediation test
mediation <- mediate(medmodel1, # Mediator model
                    outputmodel1, # Outcome model
                    boot = T, # Ask for bootstrapped confidence intervals
                    treat="emo_health_index", # Name of the x variable
                    mediator="frsbe_total" # Name of the m variable
                    )
# if you don't want bootstrap, just delete 'sims' line and set boot = F

summary(mediation)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value  
## ACME            0.174001     0.000568         0.34   0.046 *
## ADE             0.200114    -0.182999         0.63   0.300  
## Total Effect    0.374116     0.020563         0.73   0.038 *
## Prop. Mediated  0.465100    -0.113592         2.70   0.084 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 45 
## 
## 
## Simulations: 1000

plot(mediation)

There is no significant indirect effect indicating no mediation

patha <- paste(round(as.numeric(coef(medmodel1)[2]),2))
patha

## [1] "-0.69"

pathb <- paste(round(as.numeric(coef(outputmodel1)[3]),2))
pathb

## [1] "-0.25"

# indirect
direct <- paste(round(mediation$z.avg,2)) # direct
total <- paste(round((mediation$d.avg + mediation$z.avg),2)) # total

pathc <- paste(direct,"(",total,")")

data <- c(0, patha, 0,
          0, 0, 0, 
          pathb, pathc, 0) # order: top = effect of IV on mediator; bottom left = effect of mediator on DV; bottom middle = direct effect (total effect)
M<- matrix (nrow=3, ncol=3, byrow = TRUE, data=data)
plot<- plotmat (M, pos=c(1,2), 
                name= c("Self-regulation: Total","Emotional Health", "ACS Global"), # order: Mediator, IV, DV 
                box.type = "rect", box.size = 0.12, box.prop=0.5,  curve=0)

#### SPS

# Initial Model
model1 <- lm(acsg_retain ~ spstotal, df) # Y ~ X, DV predicted by IV - no mediation considered - total effect
summary(model1)

## 
## Call:
## lm(formula = acsg_retain ~ spstotal, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.859  -8.511  -0.490   7.444  31.792 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  39.5176    16.9682   2.329   0.0246 *
## spstotal      0.4564     0.2114   2.160   0.0364 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.96 on 43 degrees of freedom
## Multiple R-squared:  0.09785,    Adjusted R-squared:  0.07687 
## F-statistic: 4.664 on 1 and 43 DF,  p-value: 0.03643

# Mediation paths
medmodel1 <- lm(frsbe_total ~ spstotal, df) # M ~ X, mediator predicted by X
outputmodel1 <- lm(acsg_retain ~ spstotal + frsbe_total, df) # Y ~ X + M, DV predicted by mediator, adjusting for IV

summary(medmodel1)

## 
## Call:
## lm(formula = frsbe_total ~ spstotal, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.221 -12.092  -2.445  10.451  36.709 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 176.2946    18.4017   9.580 3.12e-12 ***
## spstotal     -0.8706     0.2292  -3.798 0.000453 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.3 on 43 degrees of freedom
## Multiple R-squared:  0.2512, Adjusted R-squared:  0.2338 
## F-statistic: 14.43 on 1 and 43 DF,  p-value: 0.0004528

summary(outputmodel1)

## 
## Call:
## lm(formula = acsg_retain ~ spstotal + frsbe_total, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.387  -9.190   0.746   7.931  30.734 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  83.9669    29.2386   2.872  0.00637 **
## spstotal      0.2369     0.2377   0.997  0.32464   
## frsbe_total  -0.2521     0.1369  -1.842  0.07251 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.53 on 42 degrees of freedom
## Multiple R-squared:  0.1653, Adjusted R-squared:  0.1255 
## F-statistic: 4.159 on 2 and 42 DF,  p-value: 0.0225

# Mediation test
mediation <- mediate(medmodel1, # Mediator model
                    outputmodel1, # Outcome model
                    boot = T, # Ask for bootstrapped confidence intervals
                    treat="spstotal", # Name of the x variable
                    mediator="frsbe_total" # Name of the m variable
                    )
# if you don't want bootstrap, just delete 'sims' line and set boot = F

summary(mediation)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value   
## ACME             0.2195       0.0469         0.49   0.008 **
## ADE              0.2369      -0.2789         0.67   0.364   
## Total Effect     0.4564       0.0261         0.83   0.036 * 
## Prop. Mediated   0.4809       0.0345         3.35   0.044 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 45 
## 
## 
## Simulations: 1000

plot(mediation)

There is a significant indirect effect and insignificant direct effect indicating full mediation

patha <- paste(round(as.numeric(coef(medmodel1)[2]),2))
patha

## [1] "-0.87"

pathb <- paste(round(as.numeric(coef(outputmodel1)[3]),2))
pathb

## [1] "-0.25"

# indirect
direct <- paste(round(mediation$z.avg,2)) # direct
total <- paste(round((mediation$d.avg + mediation$z.avg),2)) # total

pathc <- paste(direct,"(",total,")")

data <- c(0, patha, 0,
          0, 0, 0, 
          pathb, pathc, 0) # order: top = effect of IV on mediator; bottom left = effect of mediator on DV; bottom middle = direct effect (total effect)
M<- matrix (nrow=3, ncol=3, byrow = TRUE, data=data)
plot<- plotmat (M, pos=c(1,2), 
                name= c("Self-regulation: Total","Social support", "ACS Global"), # order: Mediator, IV, DV 
                box.type = "rect", box.size = 0.12, box.prop=0.5,  curve=0)

Post Hoc Mediation

In this set of mediation analysis, we look specifically at apathy as a mediator for the ACS social re-engagement outcome

# Initial Model
model1 <- lm(acss_retain ~ tbiqol_depression_tscore, df) # Y ~ X, DV predicted by IV - no mediation considered - total effect
summary(model1)

## 
## Call:
## lm(formula = acss_retain ~ tbiqol_depression_tscore, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.291 -10.478  -0.695  11.332  34.092 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              116.4189    14.9093   7.808 8.87e-10 ***
## tbiqol_depression_tscore  -0.7946     0.2735  -2.905  0.00578 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.58 on 43 degrees of freedom
## Multiple R-squared:  0.1641, Adjusted R-squared:  0.1446 
## F-statistic: 8.441 on 1 and 43 DF,  p-value: 0.005776

# Mediation paths
medmodel1 <- lm(frsbe_apathy ~ tbiqol_depression_tscore, df) # M ~ X, mediator predicted by X
outputmodel1 <- lm(acss_retain ~ tbiqol_depression_tscore + frsbe_apathy, df) # Y ~ X + M, DV predicted by mediator, adjusting for IV

summary(medmodel1)

## 
## Call:
## lm(formula = frsbe_apathy ~ tbiqol_depression_tscore, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.476  -5.933  -2.038   5.234  21.086 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)   
## (Intercept)               12.7754     6.2857   2.032  0.04831 * 
## tbiqol_depression_tscore   0.3740     0.1153   3.243  0.00229 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.412 on 43 degrees of freedom
## Multiple R-squared:  0.1965, Adjusted R-squared:  0.1779 
## F-statistic: 10.52 on 1 and 43 DF,  p-value: 0.002288

summary(outputmodel1)

## 
## Call:
## lm(formula = acss_retain ~ tbiqol_depression_tscore + frsbe_apathy, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.310 -10.881  -1.689  15.052  26.769 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              127.1919    14.7620   8.616 7.86e-11 ***
## tbiqol_depression_tscore  -0.4793     0.2886  -1.661   0.1042    
## frsbe_apathy              -0.8433     0.3421  -2.465   0.0179 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.63 on 42 degrees of freedom
## Multiple R-squared:  0.2697, Adjusted R-squared:  0.235 
## F-statistic: 7.757 on 2 and 42 DF,  p-value: 0.001359

# Mediation test
mediation <- mediate(medmodel1, # Mediator model
                    outputmodel1, # Outcome model
                    boot = T, # Ask for bootstrapped confidence intervals
                    treat="tbiqol_depression_tscore", # Name of the x variable
                    mediator="frsbe_apathy" # Name of the m variable
                    )
# if you don't want bootstrap, just delete 'sims' line and set boot = F

summary(mediation)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME            -0.3154      -0.7415        -0.05   0.014 *  
## ADE             -0.4793      -1.0666         0.18   0.142    
## Total Effect    -0.7946      -1.3003        -0.26  <2e-16 ***
## Prop. Mediated   0.3969       0.0499         1.49   0.014 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 45 
## 
## 
## Simulations: 1000

plot(mediation)

patha <- paste(round(as.numeric(coef(medmodel1)[2]),2))
patha

## [1] "0.37"

pathb <- paste(round(as.numeric(coef(outputmodel1)[3]),2))
pathb

## [1] "-0.84"

# indirect
direct <- paste(round(mediation$z.avg,2)) # direct
total <- paste(round((mediation$d.avg + mediation$z.avg),2)) # total

pathc <- paste(direct,"(",total,")")

data <- c(0, patha, 0,
          0, 0, 0, 
          pathb, pathc, 0) # order: top = effect of IV on mediator; bottom left = effect of mediator on DV; bottom middle = direct effect (total effect)
M<- matrix (nrow=3, ncol=3, byrow = TRUE, data=data)
plot<- plotmat (M, pos=c(1,2), 
                name= c("Self-regulation: Apathy","Depression", "ACS Social"), # order: Mediator, IV, DV 
                box.type = "rect", box.size = 0.12, box.prop=0.5,  curve=0)

# Initial Model
model1 <- lm(acss_retain ~ tbiqol_fatigue_tscore, df) # Y ~ X, DV predicted by IV - no mediation considered - total effect
summary(model1)

## 
## Call:
## lm(formula = acss_retain ~ tbiqol_fatigue_tscore, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.816 -11.247   0.976  14.487  39.910 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           128.2512    16.9075   7.585 1.85e-09 ***
## tbiqol_fatigue_tscore  -0.9953     0.3054  -3.260  0.00218 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.22 on 43 degrees of freedom
## Multiple R-squared:  0.1981, Adjusted R-squared:  0.1795 
## F-statistic: 10.63 on 1 and 43 DF,  p-value: 0.002185

# Mediation paths
medmodel1 <- lm(frsbe_apathy ~ tbiqol_fatigue_tscore, df) # M ~ X, mediator predicted by X
outputmodel1 <- lm(acss_retain ~ tbiqol_fatigue_tscore + frsbe_apathy, df) # Y ~ X + M, DV predicted by mediator, adjusting for IV

summary(medmodel1)

## 
## Call:
## lm(formula = frsbe_apathy ~ tbiqol_fatigue_tscore, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1646  -4.9697  -0.0723   4.2892  16.6225 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             7.5055     7.1163   1.055 0.297460    
## tbiqol_fatigue_tscore   0.4630     0.1285   3.602 0.000812 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.247 on 43 degrees of freedom
## Multiple R-squared:  0.2318, Adjusted R-squared:  0.214 
## F-statistic: 12.98 on 1 and 43 DF,  p-value: 0.0008123

summary(outputmodel1)

## 
## Call:
## lm(formula = acss_retain ~ tbiqol_fatigue_tscore + frsbe_apathy, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.425 -13.297  -1.765  13.642  39.854 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           134.0816    16.3751   8.188 3.06e-10 ***
## tbiqol_fatigue_tscore  -0.6357     0.3331  -1.908   0.0632 .  
## frsbe_apathy           -0.7768     0.3465  -2.242   0.0303 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.46 on 42 degrees of freedom
## Multiple R-squared:  0.2839, Adjusted R-squared:  0.2498 
## F-statistic: 8.324 on 2 and 42 DF,  p-value: 0.0009015

# Mediation test
mediation <- mediate(medmodel1, # Mediator model
                    outputmodel1, # Outcome model
                    boot = T, # Ask for bootstrapped confidence intervals
                    treat="tbiqol_fatigue_tscore", # Name of the x variable
                    mediator="frsbe_apathy" # Name of the m variable
                    )
# if you don't want bootstrap, just delete 'sims' line and set boot = F

summary(mediation)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value   
## ACME            -0.3597      -0.8362        -0.10   0.002 **
## ADE             -0.6357      -1.2806         0.07   0.074 . 
## Total Effect    -0.9953      -1.6721        -0.42   0.002 **
## Prop. Mediated   0.3613       0.0964         1.12   0.004 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 45 
## 
## 
## Simulations: 1000

plot(mediation)

patha <- paste(round(as.numeric(coef(medmodel1)[2]),2))
patha

## [1] "0.46"

pathb <- paste(round(as.numeric(coef(outputmodel1)[3]),2))
pathb

## [1] "-0.78"

# indirect
direct <- paste(round(mediation$z.avg,2)) # direct
total <- paste(round((mediation$d.avg + mediation$z.avg),2)) # total

pathc <- paste(direct,"(",total,")")

data <- c(0, patha, 0,
          0, 0, 0, 
          pathb, pathc, 0) # order: top = effect of IV on mediator; bottom left = effect of mediator on DV; bottom middle = direct effect (total effect)
M<- matrix (nrow=3, ncol=3, byrow = TRUE, data=data)
plot<- plotmat (M, pos=c(1,2), 
                name= c("Self-regulation: Apathy","Fatigue", "ACS Social"), # order: Mediator, IV, DV 
                box.type = "rect", box.size = 0.12, box.prop=0.5,  curve=0)

# Initial Model
model1 <- lm(acss_retain ~ tbiqol_genconcern_tscore, df) # Y ~ X, DV predicted by IV - no mediation considered - total effect
summary(model1)

## 
## Call:
## lm(formula = acss_retain ~ tbiqol_genconcern_tscore, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.441 -12.621   3.836  10.889  26.516 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                31.329     10.717   2.923 0.005503 ** 
## tbiqol_genconcern_tscore    1.163      0.286   4.067 0.000199 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.34 on 43 degrees of freedom
## Multiple R-squared:  0.2778, Adjusted R-squared:  0.2611 
## F-statistic: 16.54 on 1 and 43 DF,  p-value: 0.0001991

# Mediation paths
medmodel1 <- lm(frsbe_apathy ~ tbiqol_genconcern_tscore, df) # M ~ X, mediator predicted by X
outputmodel1 <- lm(acss_retain ~ tbiqol_genconcern_tscore + frsbe_apathy, df) # Y ~ X + M, DV predicted by mediator, adjusting for IV

summary(medmodel1)

## 
## Call:
## lm(formula = frsbe_apathy ~ tbiqol_genconcern_tscore, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.318  -5.181   0.571   4.978  14.688 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               46.0586     5.0127   9.188 1.06e-11 ***
## tbiqol_genconcern_tscore  -0.3621     0.1338  -2.707   0.0097 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.643 on 43 degrees of freedom
## Multiple R-squared:  0.1456, Adjusted R-squared:  0.1257 
## F-statistic: 7.328 on 1 and 43 DF,  p-value: 0.009701

summary(outputmodel1)

## 
## Call:
## lm(formula = acss_retain ~ tbiqol_genconcern_tscore + frsbe_apathy, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.918 -10.873  -0.563  12.551  29.280 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               65.1511    17.5308   3.716 0.000591 ***
## tbiqol_genconcern_tscore   0.8973     0.2940   3.052 0.003932 ** 
## frsbe_apathy              -0.7343     0.3098  -2.370 0.022446 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.53 on 42 degrees of freedom
## Multiple R-squared:  0.363,  Adjusted R-squared:  0.3327 
## F-statistic: 11.97 on 2 and 42 DF,  p-value: 7.697e-05

# Mediation test
mediation <- mediate(medmodel1, # Mediator model
                    outputmodel1, # Outcome model
                    boot = T, # Ask for bootstrapped confidence intervals
                    treat="tbiqol_genconcern_tscore", # Name of the x variable
                    mediator="frsbe_apathy" # Name of the m variable
                    )
# if you don't want bootstrap, just delete 'sims' line and set boot = F

summary(mediation)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME             0.2659       0.0253         0.59   0.030 *  
## ADE              0.8973       0.3336         1.45   0.008 ** 
## Total Effect     1.1633       0.5921         1.69  <2e-16 ***
## Prop. Mediated   0.2286       0.0188         0.54   0.030 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 45 
## 
## 
## Simulations: 1000

plot(mediation)

patha <- paste(round(as.numeric(coef(medmodel1)[2]),2))
patha

## [1] "-0.36"

pathb <- paste(round(as.numeric(coef(outputmodel1)[3]),2))
pathb

## [1] "-0.73"

# indirect
direct <- paste(round(mediation$z.avg,2)) # direct
total <- paste(round((mediation$d.avg + mediation$z.avg),2)) # total

pathc <- paste(direct,"(",total,")")

data <- c(0, patha, 0,
          0, 0, 0, 
          pathb, pathc, 0) # order: top = effect of IV on mediator; bottom left = effect of mediator on DV; bottom middle = direct effect (total effect)
M<- matrix (nrow=3, ncol=3, byrow = TRUE, data=data)
plot<- plotmat (M, pos=c(1,2), 
                name= c("Self-regulation: Apathy","General cognition", "ACS Social"), # order: Mediator, IV, DV 
                box.type = "rect", box.size = 0.12, box.prop=0.5,  curve=0)

RQ2: Time Since Injury

2. How does time since injury influence resiliency?

1. Is time after injury associated with an individual’s ability to re-engage in meaningful activities?
2. To what extent do the relationships between protective factors and self-regulatory processes with resiliency-related outcomes change with time since sustaining TBI? Hypothesis: Time since injury will be associated with resiliency, leading to different resiliency-related outcomes

To answer these questions, first look at descriptive statistics, then regression model with time since injury included, lastly, investigate what, if any, role time since injury has on protective factors and self-regulation

Descriptives

selected_vars_acs <- c("time_injury")

summary_table <- do.call(rbind, lapply(df[selected_vars_acs], function(x) round(summary(x), 2)))

print(summary_table)

##             Min. 1st Qu. Median Mean 3rd Qu. Max.
## time_injury    1     3.5      7 8.91      11   30

## [1] 8.911111

## [1] 7.07933

Correlation ACS

# Define variables
timesinceinjury_variables <- c("TimeSinceInjury", "Global", "IADL", "Leisure", "Fitness", "Social", "acsg_before", "acsg_after", "acsi_before", "acsi_after","acsl_before", "acsl_after","acsf_before", "acsf_after","acss_before", "acss_after")

# Ensure selected variables are in the dataframe
timesinceinjury_variables <- intersect(timesinceinjury_variables, colnames(df2))

# Extract relevant data
timesinceinjury_df2 <- df2[, timesinceinjury_variables]

# Calculate correlation matrix
cor_matrix <- rcorr(as.matrix(timesinceinjury_df2), type = "spearman")$r
cor_matrix[upper.tri(cor_matrix)] <- NA
p_matrix <- rcorr(as.matrix(timesinceinjury_df2), type = "spearman")$P
p_matrix[is.na(p_matrix)] <- .0000001
p_matrix[upper.tri(p_matrix)] <- NA

# Melt the correlation matrix for ggplot
melted_cor <- melt(cor_matrix, na.rm = TRUE)
melted_p <- melt(p_matrix, na.rm = TRUE)

melted_cor$p <- melted_p$value
melted_cor$psig <- ""
melted_cor$psig[melted_cor$p < .05] <- "*"
melted_cor$psig[melted_cor$p < .01] <- "**"
melted_cor$psig[melted_cor$p < .001] <- "***"

# Create a heatmap using ggplot2
ggplot(melted_cor, aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(value,2)), vjust = 1, size = 6) +  # Adjust size here
  geom_text(aes(label = psig), vjust = .25, size = 6) +  # Adjust size here
  scale_fill_gradient2(low = "purple", mid = "white", high = "orange", 
                       midpoint = 0, limit = c(-1,1), space = "Lab",
                       name = "Correlation") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 20),  # Adjust size here
        axis.text.y = element_text(angle = 45, hjust = 1, size = 20),  # Adjust size here
        axis.title = element_text(size = 14),
        axis.ticks = element_line(size = 1)) +
  labs(caption = "<.05 = *, <.01 = **, <.001 = ***") +
  xlab("") +
  ylab("") +
  ggtitle("Correlation Matrix Personal Protective Factors with ACS Variables") +
  guides(fill = guide_colorbar(title.position = "right"))

Regression Analysis

Global ACS

First, looking just at the relationship between time since injury and re-engagement while controlling for age.

model_timeinjury <- lm(acsg_retain~age_current+time_injury, data=df)
summary(model_timeinjury)

## 
## Call:
## lm(formula = acsg_retain ~ age_current + time_injury, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.816 -10.974  -1.151   6.725  40.839 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  89.6585     8.1704  10.974 6.52e-14 ***
## age_current  -0.4016     0.1648  -2.437   0.0191 *  
## time_injury   0.5741     0.3387   1.695   0.0974 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.65 on 42 degrees of freedom
## Multiple R-squared:  0.1528, Adjusted R-squared:  0.1124 
## F-statistic: 3.787 on 2 and 42 DF,  p-value: 0.03076

Controlling for age, there is no significant relationship between time since injury and re-engagement

Social ACS3

model_timeinjury_soc <- lm(acss_retain~age_current+time_injury, data=df)
summary(model_timeinjury_soc)

## 
## Call:
## lm(formula = acss_retain ~ age_current + time_injury, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.840 -11.829   3.233  13.604  26.265 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  87.7492     9.1450   9.595  3.8e-12 ***
## age_current  -0.4611     0.1845  -2.499   0.0164 *  
## time_injury   0.8767     0.3791   2.313   0.0257 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.51 on 42 degrees of freedom
## Multiple R-squared:  0.1898, Adjusted R-squared:  0.1512 
## F-statistic: 4.919 on 2 and 42 DF,  p-value: 0.01204

IADL ACS3

model_timeinjury_iadl <- lm(acsi_retain~age_current+time_injury, data=df)
summary(model_timeinjury_iadl)

## 
## Call:
## lm(formula = acsi_retain ~ age_current + time_injury, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.051 -14.366  -1.745  15.838  31.228 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  89.9029     9.3742   9.590 3.85e-12 ***
## age_current  -0.3031     0.1891  -1.603    0.116    
## time_injury   0.5357     0.3886   1.379    0.175    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.95 on 42 degrees of freedom
## Multiple R-squared:  0.08296,    Adjusted R-squared:  0.03929 
## F-statistic:   1.9 on 2 and 42 DF,  p-value: 0.1622

Leisure ACS3

model_timeinjury_leisure <- lm(acsl_retain~age_current+time_injury, data=df)
summary(model_timeinjury_leisure)

## 
## Call:
## lm(formula = acsl_retain ~ age_current + time_injury, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -28.48 -12.54   0.29  12.55  53.19 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  97.0092     9.2193  10.522  2.4e-13 ***
## age_current  -0.4573     0.1860  -2.459   0.0181 *  
## time_injury   0.5682     0.3821   1.487   0.1445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.65 on 42 degrees of freedom
## Multiple R-squared:  0.1459, Adjusted R-squared:  0.1053 
## F-statistic: 3.588 on 2 and 42 DF,  p-value: 0.03641

Fitness ACS3

model_timeinjury_fitness <- lm(acsf_retain~age_current+time_injury, data=df)
summary(model_timeinjury_fitness)

## 
## Call:
## lm(formula = acsf_retain ~ age_current + time_injury, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.880 -19.191  -5.599   5.963 129.732 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  80.0033    17.2377   4.641 3.38e-05 ***
## age_current  -0.4660     0.3477  -1.340    0.187    
## time_injury   0.9567     0.7145   1.339    0.188    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.01 on 42 degrees of freedom
## Multiple R-squared:  0.06756,    Adjusted R-squared:  0.02316 
## F-statistic: 1.522 on 2 and 42 DF,  p-value: 0.2301

Hierarchical Regression: Step 4

We add time since injury into the previous hierarchical regression model to see if it is a significant predictor of recovery

step4<- lm(acsg_retain~age_current+phys_health_index+tbiqol_genconcern_tscore+emo_health_index+ spstotal+frsbe_total+time_injury, data=df)
summary(step4)

## 
## Call:
## lm(formula = acsg_retain ~ age_current + phys_health_index + 
##     tbiqol_genconcern_tscore + emo_health_index + spstotal + 
##     frsbe_total + time_injury, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.870  -7.846  -1.752   6.897  32.012 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)               48.8467    36.8221   1.327   0.1928  
## age_current               -0.1402     0.1651  -0.849   0.4012  
## phys_health_index          0.1009     0.2204   0.458   0.6499  
## tbiqol_genconcern_tscore   0.7836     0.3859   2.031   0.0495 *
## emo_health_index          -0.1265     0.2172  -0.582   0.5638  
## spstotal                   0.1809     0.2447   0.739   0.4645  
## frsbe_total               -0.1167     0.1381  -0.845   0.4036  
## time_injury                0.6899     0.3094   2.230   0.0319 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.93 on 37 degrees of freedom
## Multiple R-squared:  0.408,  Adjusted R-squared:  0.2961 
## F-statistic: 3.643 on 7 and 37 DF,  p-value: 0.00436

#Nested Model Comparison
anova(step3, step4)

## Analysis of Variance Table
## 
## Model 1: acsg_retain ~ age_current + phys_health_index + tbiqol_genconcern_tscore + 
##     emo_health_index + spstotal + frsbe_total
## Model 2: acsg_retain ~ age_current + phys_health_index + tbiqol_genconcern_tscore + 
##     emo_health_index + spstotal + frsbe_total + time_injury
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     38 8148.8                              
## 2     37 7183.5  1    965.25 4.9717 0.03192 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#change in R-squared
summary(step4)$r.squared - summary(step3)$r.squared

## [1] 0.0795411

POST HOC: Groups

Because there was no significant relationship between time since injury and outcomes BUT it was a significant predictor in the model, I seperated time since recovery into three groups to examine relationships further. Early = 6mo-3years Mid = 3 to 10 years Later = >10 years

1= everything <=3 and 3 is everything >=10

demo_table_time <- df_time3 %>%
  subset(., select = c(time_injury_ex, age_current, work_current, substance, severity, acsg_retain, acss_retain, acsi_retain, acsl_retain, acsf_retain)) %>%
  tbl_summary(
    missing = "no",
    by = time_injury_ex,
    type = list(
      c(age_current, substance, acsg_retain, acss_retain, acsi_retain, acsl_retain, acsf_retain) ~ "continuous",
      c(work_current, severity) ~ "categorical"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} ({p}%)"),
    label = list(
      age_current ~ "Age (years)",
      time_injury_ex ~ "Early/Late Recovery",
      work_current ~ "Employment status",
      substance ~ "Substance use score",
      severity ~ "Severity of Injury",
      acsg_retain ~ "Global ACS",
      acss_retain ~ "Social ACS",
      acsi_retain ~ "IADL ACS",
      acsl_retain ~ "Leisure ACS",
      acsf_retain ~ "Fitness ACS"
    )
  ) 


print(demo_table_time)   

Characteristic	Early, N = 111	Mid, N = 141	Later, N = 201
Age (years)	41 (13)	50 (17)	49 (13)
Employment status
0	4 (36%)	11 (79%)	11 (55%)
1	7 (64%)	3 (21%)	9 (45%)
Substance use score	1.82 (1.66)	3.21 (3.53)	2.95 (3.17)
Severity of Injury
2	5 (45%)	7 (50%)	7 (35%)
3	6 (55%)	7 (50%)	13 (65%)
Global ACS	71 (14)	73 (18)	80 (16)
Social ACS	70 (17)	69 (14)	79 (22)
IADL ACS	74 (19)	79 (20)	85 (16)
Leisure ACS	74 (17)	79 (21)	85 (17)
Fitness ACS	61 (21)	61 (31)	73 (40)
1 Mean (SD); n (%)

level_counts3 <- table(df_time3$time_injury_ex)

# Displaying the counts
print(level_counts3)

## 
## Early   Mid Later 
##    11    14    20

We see the counts of # participants in each group

Global

Global ACS3 scores (ie, global re-engagement scores)

group_by(df_time3, time_injury_ex) %>%
  summarise(
    count = n(),
    mean = mean(acsg_retain, na.rm = TRUE),
    sd = sd(acsg_retain, na.rm = TRUE)
  )

## # A tibble: 3 × 4
##   time_injury_ex count  mean    sd
##   <fct>          <int> <dbl> <dbl>
## 1 Early             11  70.9  13.9
## 2 Mid               14  73.4  18.3
## 3 Later             20  80.2  16.4

library("ggpubr")
ggboxplot(df_time3, x = "time_injury_ex", y = "acsg_retain",
          color = "time_injury_ex", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Early", "Mid", "Later"),
          ylab = "Global Engagement", xlab = "Time Since Injury")

library("ggpubr")

# Set global font to Times New Roman
theme_set(theme_gray(base_family = "Times New Roman"))

# Create the plot
plot <- ggline(df_time3, x = "time_injury_ex", y = "acsg_retain",
               add = c("mean_se", "jitter"),
               order = c("Early", "Mid", "Late")) +
        labs(x = "Time Since Injury", y = "Global Engagement", 
             face = "bold", family = "Times New Roman")

# Save the plot as a PNG file with DPI=300
ggsave("plot.png", plot, dpi = 300)

library("ggpubr")
ggline(df_time3, x = "time_injury_ex", y = "acsg_retain",
       add = c("mean_se", "jitter"),
       order = c("Early", "Mid", "Later"),
       ylab = "Global Engagement", xlab = "Time Since Injury")

res.aov <- aov(acsg_retain~time_injury_ex, data=df_time3)
summary(res.aov)

##                Df Sum Sq Mean Sq F value Pr(>F)
## time_injury_ex  2    734   366.9   1.352   0.27
## Residuals      42  11401   271.5

model_timeinjury_ex <- lm(acsg_retain~age_current+time_injury_ex, data=df_time3)
summary(model_timeinjury_ex)

## 
## Call:
## lm(formula = acsg_retain ~ age_current + time_injury_ex, data = df_time3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.122 -11.225  -0.144   9.222  39.345 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          88.0833     8.2434  10.685 2.02e-13 ***
## age_current          -0.4207     0.1663  -2.530   0.0154 *  
## time_injury_exMid     6.3414     6.4355   0.985   0.3302    
## time_injury_exLater  12.6493     5.9712   2.118   0.0403 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.51 on 41 degrees of freedom
## Multiple R-squared:  0.1873, Adjusted R-squared:  0.1279 
## F-statistic:  3.15 on 3 and 41 DF,  p-value: 0.03502

Controlling for age, we see a significant relationship between time since injury and global re engagement- specifically between early and later recovery

Social

social ACS3 scores (ie, social re-engagement scores)

group_by(df_time3, time_injury_ex) %>%
  summarise(
    count = n(),
    mean = mean(acss_retain, na.rm = TRUE),
    sd = sd(acss_retain, na.rm = TRUE)
  )

## # A tibble: 3 × 4
##   time_injury_ex count  mean    sd
##   <fct>          <int> <dbl> <dbl>
## 1 Early             11  69.8  17.3
## 2 Mid               14  69.4  14.2
## 3 Later             20  79    22.2

library("ggpubr")
ggboxplot(df_time3, x = "time_injury_ex", y = "acss_retain",
          color = "time_injury_ex", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Early", "Mid", "Later"),
          ylab = "Social Engagement", xlab = "Time Since Injury")

library("ggpubr")
ggline(df_time3, x = "time_injury_ex", y = "acss_retain",
       add = c("mean_se", "jitter"),
       order = c("Early", "Mid", "Later"),
       ylab = "Social Engagement", xlab = "Time Since Injury")

res.aov <- aov(acss_retain~time_injury_ex, data=df_time3)
summary(res.aov)

##                Df Sum Sq Mean Sq F value Pr(>F)
## time_injury_ex  2    983   491.4   1.384  0.262
## Residuals      42  14915   355.1

model_timeinjury_ex <- lm(acss_retain~age_current+time_injury_ex, data=df_time3)
summary(model_timeinjury_ex)

## 
## Call:
## lm(formula = acss_retain ~ age_current + time_injury_ex, data = df_time3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.838 -12.003   4.362  14.446  27.731 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          87.8939     9.5406   9.213 1.54e-11 ***
## age_current          -0.4428     0.1925  -2.301   0.0266 *  
## time_injury_exMid     3.7081     7.4482   0.498   0.6213    
## time_injury_exLater  12.7164     6.9109   1.840   0.0730 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.95 on 41 degrees of freedom
## Multiple R-squared:  0.1691, Adjusted R-squared:  0.1083 
## F-statistic: 2.781 on 3 and 41 DF,  p-value: 0.05301

IADL

IADL ACS3 scores (ie, IADL re-engagement scores)

group_by(df_time3, time_injury_ex) %>%
  summarise(
    count = n(),
    mean = mean(acsi_retain, na.rm = TRUE),
    sd = sd(acsi_retain, na.rm = TRUE)
  )

## # A tibble: 3 × 4
##   time_injury_ex count  mean    sd
##   <fct>          <int> <dbl> <dbl>
## 1 Early             11  73.5  19.1
## 2 Mid               14  79.1  20.2
## 3 Later             20  85.0  16.0

library("ggpubr")
ggboxplot(df_time3, x = "time_injury_ex", y = "acsi_retain",
          color = "time_injury_ex", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Early", "Mid", "Later"),
          ylab = "IADL Engagement", xlab = "Time Since Injury")

library("ggpubr")
ggline(df_time3, x = "time_injury_ex", y = "acsi_retain",
       add = c("mean_se", "jitter"),
       order = c("Early", "Mid", "Later"),
       ylab = "IADL Engagement", xlab = "Time Since Injury")

ANOVA to see if there is a difference between the 3 time groups

res.aov <- aov(acsi_retain~time_injury_ex, data=df_time3)
summary(res.aov)

##                Df Sum Sq Mean Sq F value Pr(>F)
## time_injury_ex  2    953   476.5    1.45  0.246
## Residuals      42  13805   328.7

This is controlling for age to see if there is a difference between the three groups

model_timeinjury_ex <- lm(acsi_retain~age_current+time_injury_ex, data=df_time3)
summary(model_timeinjury_ex)

## 
## Call:
## lm(formula = acsi_retain ~ age_current + time_injury_ex, data = df_time3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -28.0758 -14.6303  -0.8689  15.3251  29.2851 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          87.4550     9.3898   9.314 1.13e-11 ***
## age_current          -0.3408     0.1894  -1.799   0.0794 .  
## time_injury_exMid     8.7506     7.3305   1.194   0.2394    
## time_injury_exLater  14.1245     6.8017   2.077   0.0441 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.67 on 41 degrees of freedom
## Multiple R-squared:  0.133,  Adjusted R-squared:  0.06955 
## F-statistic: 2.096 on 3 and 41 DF,  p-value: 0.1155

We see a significant difference between early and late groups with IADl engagement when controlling for age

Leisure

Leisure ACS3 scores (ie, leisure re-engagement scores)

group_by(df_time3, time_injury_ex) %>%
  summarise(
    count = n(),
    mean = mean(acsl_retain, na.rm = TRUE),
    sd = sd(acsl_retain, na.rm = TRUE)
  )

## # A tibble: 3 × 4
##   time_injury_ex count  mean    sd
##   <fct>          <int> <dbl> <dbl>
## 1 Early             11  74.3  17.0
## 2 Mid               14  78.8  21.0
## 3 Later             20  85.0  17.5

library("ggpubr")
ggboxplot(df_time3, x = "time_injury_ex", y = "acsl_retain",
          color = "time_injury_ex", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Early", "Mid", "Later"),
          ylab = "Leisure Engagement", xlab = "Time Since Injury")

library("ggpubr")
ggline(df_time3, x = "time_injury_ex", y = "acsl_retain",
       add = c("mean_se", "jitter"),
       order = c("Early", "Mid", "Later"),
       ylab = "Leisure Engagement", xlab = "Time Since Injury")

res.aov <- aov(acsl_retain~time_injury_ex, data=df_time3)
summary(res.aov)

##                Df Sum Sq Mean Sq F value Pr(>F)
## time_injury_ex  2    882   440.9   1.282  0.288
## Residuals      42  14445   343.9

model_timeinjury_ex <- lm(acsl_retain~age_current+time_injury_ex, data=df_time3)
summary(model_timeinjury_ex)

## 
## Call:
## lm(formula = acsl_retain ~ age_current + time_injury_ex, data = df_time3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.534 -15.551   0.748  10.588  51.109 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          94.5042     9.2096  10.261 6.84e-13 ***
## age_current          -0.4956     0.1858  -2.667   0.0109 *  
## time_injury_exMid     9.0993     7.1899   1.266   0.2128    
## time_injury_exLater  14.7334     6.6711   2.209   0.0329 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.33 on 41 degrees of freedom
## Multiple R-squared:  0.1969, Adjusted R-squared:  0.1381 
## F-statistic: 3.351 on 3 and 41 DF,  p-value: 0.02803

We see a significant difference between early and late groups with Leisure engagement when controlling for age

Fitness

Fitness ACS3 scores (ie, fitness re-engagement scores)

group_by(df_time3, time_injury_ex) %>%
  summarise(
    count = n(),
    mean = mean(acsf_retain, na.rm = TRUE),
    sd = sd(acsf_retain, na.rm = TRUE)
  )

## # A tibble: 3 × 4
##   time_injury_ex count  mean    sd
##   <fct>          <int> <dbl> <dbl>
## 1 Early             11  60.7  20.6
## 2 Mid               14  61.4  30.6
## 3 Later             20  73.3  40.4

library("ggpubr")
ggboxplot(df_time3, x = "time_injury_ex", y = "acsg_retain",
          color = "time_injury_ex", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Early", "Mid", "Later"),
          ylab = "Fitness Engagement", xlab = "Time Since Injury")

library("ggpubr")
ggline(df_time3, x = "time_injury_ex", y = "acsf_retain",
       add = c("mean_se", "jitter"),
       order = c("Early", "Mid", "Later"),
       ylab = "Fitness Engagement", xlab = "Time Since Injury")

res.aov <- aov(acsf_retain~time_injury_ex, data=df_time3)
summary(res.aov)

##                Df Sum Sq Mean Sq F value Pr(>F)
## time_injury_ex  2   1662   830.8   0.736  0.485
## Residuals      42  47418  1129.0

model_timeinjury_ex <- lm(acsf_retain~age_current+time_injury_ex, data=df_time3)
summary(model_timeinjury_ex)

## 
## Call:
## lm(formula = acsf_retain ~ age_current + time_injury_ex, data = df_time3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -55.517 -18.751  -6.624   9.536 133.653 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          79.3987    17.7255   4.479 5.88e-05 ***
## age_current          -0.4574     0.3576  -1.279    0.208    
## time_injury_exMid     4.8626    13.8381   0.351    0.727    
## time_injury_exLater  16.2238    12.8398   1.264    0.214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.35 on 41 degrees of freedom
## Multiple R-squared:  0.07093,    Adjusted R-squared:  0.002949 
## F-statistic: 1.043 on 3 and 41 DF,  p-value: 0.3836